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ABSTRACT 

A  mixed  variational  statement  and  corresponding  finite  element 
model  are  developed  for  an  arbitrary  plane -body  undergoing  large 
deformations  (l.e.  large  displacements,  large  rotations  and  small 
strains)  under  external  loads  using  the  updated  Lagranglan 
formulation.  The  mixed  finite  element  formulation  allows  the  nodal 
displacements  and  stresses  to  be  approximated  Independently. 

Two  algorithms  are  discussed  for  the  analysis  of  a  thin,  uniformly 
loaded  plate  with  a  circular  hole  In  contact  with  a  pin.  The  different 
algorithms  consider  the  separate  cases  of  a  rigid  pin  and  a  flexible 
pin,  and  use  different  methods  to  account  for  the  computational 
difficulties  that  arise  from  the  unknown  contact  area  and  the  presence 
of  friction  between  the  pin  and  the  plate.  A  number  of  different 
contact  problems  are  solved  using  these  two  techniques. 

A  hybrid  technique  Is  presented  that  combines  the  numerical 
technique  of  the  finite  element  method  with  the  experimental  technique 
of  moire  Interferometry.  The  displacements  at  the  pin-hole  Interface 
are  measured  from  physical  experiments  and  are  then  used  as  prescribed 
boundary  conditions  In  the  finite  element  analysis  of  the  modeled 
problem.  Results  of  this  algorithm  are  compared  with  solutions  obtained 
from  strictly  computational  algorithms  that  are  independent  of  the 
experimental  data.  The  agreement  Is  found  to  be  very  good. 
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1.  INTRODUCTION 


1.1  Motivation 

The  solution  of  contact  problems  generally  Involves  the 
determination  of  the  states  of  displacement,  strain,  and  stress  acting 
within  the  domains  of  two  or  more  separate  bodies  pressing  against  one 
another  under  external  action.  Unlike  most  problems  in  solid  and 
structural  mechanics,  where  the  critical  regions  of  interest  are  usually 
far  removed  from  the  point  of  application  of  the  loads,  the  domain  of 
interest  In  contact  problems  Is  the  area  at  or  near  the  region  of 
contact,  which  is  typically  the  region  of  load  transfer.  Contact 
stresses  are  typically  among  the  largest  stresses  found  within  the 
contacting  bodies. 

The  accurate  numerical  simulation  of  the  response  of  two  elastic 
bodies  in  contact  with  one  another  under  external  load  remains  as  one  of 
the  most  challenging  problems  of  computational  solid  mechanics  due  to 
several  Inherent  complications.  The  region  or  area  of  contact  between 
the  two  bodies  is  changing  continuously  during  the  loading,  and  is 
generally  not  known  a-priori  as  a  function  of  the  applied  loading.  In 
addition,  the  presence  of  friction  between  the  two  bodies  creates 
varying  regions  of  stick  and  slip.  Hence  this  type  of  problem  is  highly 
nonlinear,  and  although  some  exact  elasticity  solutions  exist,  these  are 
typically  for  problems  with  simple  geometrical  shapes  and  frictional 
conditions.  Numerical  techniques  must  therefore  be  used  to  solve 
problems  involving  more  complicated  geometries  and  contact  conditions. 

The  finite  element  method  has  long  been  established  as  a  versatile 
and  powerful  tool  of  analysis  for  solid  and  structural  mechanics 
problems  and  has  recently  been  applied  to  numerous  studies  of  elastic 


contact  problems.  Mixed  finite  element  models.  In  which  Independent 
approximations  of  displacements  and  stresses  are  Introduced,  have  also 
been  developed  and  applied  to  problems  In  solid  mechanics,  but  to  a  much 
lesser  extent.  Mixed  formulations  would  appear  to  be  especially 
applicable  to  contact  problems  since  the  stresses  at  the  contact 
boundary  are  computed  as  part  of  the  solution  rather  than  part  of  the 
post-computation  as  practiced  in  displacement  formulations.  These 
boundary  stresses  are  of  extreme  importance  in  computing  the  regions  of 
stick  and  slip  and  are  useful  in  other  steps  of  the  analysis  as  well. 

Among  the  many  different  types  of  contact  problems  found  in 
engineering,  the  problem  of  a  thin,  pin-loaded  plate  has  recently 
received  considerable  attention,  particularly  due  to  the  increasing  use 
of  composite  materials  in  modern  structural  applications  [1-10]. 
Elasticity  solutions  and  finite  element  approximations  have  dominated 
the  majority  of  these  analyses,  most  of  which  have  ignored  the  effects 
of  the  pin  by  assuming  that  it  is  rigid  and  hence  have  only  modeled  and 
subsequently  analyzed  the  domain  of  the  plate.  A  number  of  other 
simplifying  assumptions  were  Invoked  in  many  of  these  analyses, 
including  the  assumptions  of  a  cosinusoidal  radial  traction  on  the  hole 
boundary  as  well  as  a  constant  coefficient  of  friction  acting  between 
the  pin  and  the  plate.  In  addition,  various  assumptions  were  made  to 
approximate  the  behavior  of  the  contacting  bodies  at  the  contact  region, 
since  this  information  is  required  in  order  to  obtain  a  solution  to  the 
problem,  and  is  In  general  not  known  as  a  function  of  any  of  the 
parameters  of  the  problem. 

A  review  of  the  pertinent  literature  is  presented  in  the  following 
section  to  assess  the  gains  made  in  the  past  and  to  indicate  the 


direction  of  recent  research  efforts  In  the  analysis  of  elastic  contact 
problems.  Since  most  contact  algorithms  are  strictly  numerical,  a  brief 
review  of  geometrically  nonlinear  analysis  Is  Included,  particularly 
highlighting  the  lack  of  studies  regarding  applications  of  mixed  methods 
to  nonlinear  problems  in  plane  elasticity  in  general,  and  contact 
problems  In  particular. 

1.2  Literature  Review 

The  displacement  finite  element  model,  based  on  the  principle  of 
minimum  total  potential  energy,  has  been  applied  to  many  problems  in 
structural  engineering  and  solid  mechanics  since  the  early  1960's.  The 
mathematical  and  numerical  properties  of  the  displacement  finite  element 
method  have  been  firmly  established  and  numerous  publications  allude  to 
these  advances.  The  textbooks  by  Reddy  [11],  Zienkiewicz  [12],  and 
Bathe  [13]  contain  a  host  of  references  regarding  both  the  historical 
development  of  and  recent  advances  in  most  branches  of  finite  element 
analysis. 

Although  the  mathematical  properties  of  mixed  finite  element 
approximations  have  seen  extensive  development  in  recent  years  [14-16], 
the  applications  of  mixed  models  to  actual  physical  problems  have  been 
relatively  few.  In  linear,  two-dimensional  elasticity  problems,  the 
mixed  finite  element  equations  can  be  developed  using  the  Hellinger- 
Relssner  variational  principle  [17].  Dunham  and  Pister  [18]  were  among 
the  first  to  use  this  principle  to  Introduce  mixed  finite  element 
approximations  and  present  numerical  results  for  plane  elasticity 
problems.  They  obtained  displacements  and  stresses  that  were  superior 
to  those  obtained  using  an  equivalent  displacement  model.  Pitkaranta 
and  Sternberg  [19]  analyzed  several  mixed  finite  element  methods  for  the 


plane  elasticity  equations  on  a  rectangular  domain,  while  Mlrza  and 
Olsen  [20]  presented  a  more  thorough  study  regarding  the  convergence  and 
performance  of  the  mixed  finite  element  method  for  linear  plane 
elasticity  applications. 

Contact  problems  are  generally  regarded  as  being  highly  nonlinear 
for  two  major  reasons.  First,  the  body  or  bodies  under  analysis  might 
undergo  large  deformations,  and  hence  must  be  modeled  such  that  the  new 
geometry  and  the  most  recent  state  of  stress  In  each  body  are  accounted 
for  In  an  appropriate  manner  In  the  subsequent  loading.  Second,  the 
boundary  conditions  change  continuously  as  a  result  of  the  changes  in 
the  contact  region  with  increasing  load.  The  development  of  the 
governing  finite  element  models  for  the  nonlinear  response  of  a  solid 
body  under  external  load  has  been  developed  by  numerous  investigators 
121-32] ,  but  none  of  these  have  addressed  mixed  models  for  plane 
elasticity  problems.  Horrlgmoe  and  Bergen  [33]  presented  an  incremental 
mixed  variational  principle  and  corresponding  finite  element  model  for 
solid  bodies,  but  gave  no  numerical  examples  or  comments  on 
Implementation. 

Contact  problems  have  challenged  mathematicians  and  engineers  for 
over  a  century.  In  1881,  Hertz  obtained  a  solution  for  the  problem  of 
two  elastic  cylinders  In  contact  with  one  another  [34].  Numerous 
attempts  have  since  been  made  to  accurately  model  the  physics  of  contact 
for  more  complicated  problems  using,  among  other  methods,  finite  element 
techniques.  Due  to  the  inherent  nonlinearity  of  the  contact  problem, 
all  of  the  most  successful  of  these  algorithms  contain  several  iterative 
procedures  to  account  for  the  varying  regions  of  contact  and  of  stick 


Nearly  all  investigations  of  contact  problems  using  the  finite 
element  method  have  employed  the  more  conventional  displacement 
formulations.  Chan  and  Tuba  (35]  used  conventional  elements  and 
Iterative  procedures  to  solve  several  different  contact  problems.  Bathe 
and  Chaudhary  [36]  Imposed  the  contact  conditions  by  constructing  the 
total  potential  of  the  nodal  contact  forces  for  an  element  in  contact 
and  adding  this  term  to  the  original  potential  energy  functional. 

Campos,  Oden  and  Klkuchl  (37]  solved  discretized  contact  problems  using 
prescribed  normal  boundary  tractions  and  nonlinear  inequalities. 
Francavllla  and  Zienklewlcz  [38]  used  the  flexibility  matrices  of  two 
elastic  bodies  In  frictionless  contact  along  with  iterative  procedures 
to  check  for  penetration.  This  technique  was  later  modified  by  Sachdeva 
and  Ramakrl shnan  [39]  to  include  the  effects  of  friction.  Marks  and  his 
colleagues  [40-41]  presented  several  solution  techniques  for  contact 
problems  using  the  conjugate  gradient  technique  integrated  with  the 
finite  element  method  for  frictionless  contact  problems.  Okamoto  and 
Nakazawa  [42]  presented  a  technique  which  used  three-dimensional 
elements  and  used  the  magnitude  of  load  causing  a  change  in  the  contact 
status  of  one  node  as  a  load  step.  Fredrickson  [43]  used  iterative 
techniques  to  account  for  the  contact  conditions  and  also  used  a 
superelement  technique  to  reduce  the  number  of  degrees  of  freedom. 

Past  studies  of  contact  problems  using  mixed  finite  element  models 
are  not  nearly  as  numerous  as  those  using  displacement  models. 

Hasllnger  and  Hlavacek  [44]  presented  a  mixed  formulation  for  the 
Signorini  problem  with  prescribed  normal  contact  forces  but  gave  no 
numerical  results.  Tseng  and  Olsen  [45]  applied  the  mixed  finite 


element  method  to  several  plane  elasticity  contact  problems  using  the 
equations  of  linear  elasticity  along  with  an  iterative  scheme. 

The  increasing  use  of  composite  materials  in  structural 
applications  has  generated  a  great  deal  of  interest  in  the  problem 
Involving  a  bolt  or  pin  in  contact  with  an  elastic  plate.  Although  a 
variety  of  techniques  have  been  used  to  analyze  this  problem,  most 
studies  have  used  either  elasticity  or  finite  element  solutions  to 
obtain  the  states  of  displacement  and  stress  in  the  plate  [1-101.  Most 
analyses  have  neglected  the  elasticity  of  the  pin  by  assuming  that  it  is 
perfectly  rigid  [4-9],  while  others  have  assumed  a  cosinusoidal  radial 
stress  acting  between  the  pin  and  the  plate  at  the  points  of  contact 
[1,2,10].  The  validity  of  this  latter  assumption  has  been  demonstrated 
experimentally  for  Isotropic  plates  [3],  but  was  recently  shown  to  be 
incorrect  for  orthotropic  plates  by  Hyer  and  Klang  [46].  Of  the 
numerical  studies  of  the  pin-loaded  plate  problem,  the  only  analysis 
that  accounted  for  the  combined  effects  of  nonlinearity,  actual  boundary 
loading,  friction,  and  orthotropy  of  the  plate  was  developed  by 
Wilkinson,  Rowlands,  and  Cook  [47].  They  presented  a  simple  iterative 
technique  to  compute  the  stresses  around  the  hole  of  a  pin-loaded 
orthotropic  plate.  This  method  was  later  modified  by  Rahman  et  al. 

[48].  The  results  of  these  two  techniques  were  also  verified 
experimentally  by  Wilkinson  and  Rowlands  [49].  A  review  of  other 
pertinent  methods  in  the  displacement  finite  element  analysis  of  pin 
joints  is  given  by  Rao  [50]. 

1.3  The  Present  Study 

Despite  the  gains  made  in  recent  years  in  the  understanding  of 
contact  phenomena  between  solid  bodies,  the  extreme  complexity  of  such 
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problems  has  thwarted  the  development  of  a  general  and  effective 
computational  method  of  analysis.  In  particular,  mixed  finite  element 
models,  which  contain  stresses  as  nodal  variables  and  hence  would  seem 
to  be  particularly  suited  for  the  analysis  of  contact  problems,  have 
seen  extremely  limited  use,  particularly  in  geometrically  nonlinear 
analysis. 

Furthermore,  all  numerical  algorithms,  regardless  of  the 
formulation  used,  are  based  on  several  key  assumptions  and 
approximations  when  performing  an  analysis  on  two  or  more  bodies  in 
contact.  Various  important  quantities,  such  as  the  static  and  dynamic 
coefficients  of  friction,  are  assumed  to  be  constant  as  a  function  of 
position  and  load.  Hence  it  is  difficult  to  isolate  the  effects  of  the 
different  assumptions  made  in  most  computational  schemes  that  have  been 
developed  in  the  past. 

The  present  study  will  address  each  of  these  major  difficulties  or 
limitations  and  will  primarily  be  directed  toward  the  development  of  a 
computational  technique  for  the  analysis  of  two-dimensional  contact 
problems.  Different  strategies  developed  here  will  use  a  mixed  updated 
Lagrangian  formulation  as  a  basis  for  the  analysis.  This  formulation 
has  seen  little  or  no  use  as  a  tool  of  analysis  for  elasticity  problems 
in  the  past;  however,  it  would  be  well  suited  for  the  analysis  of 
contact  problems  since  the  displacements  and  stresses  are  approximated 
independently  and  each  of  the  respective  components  appear  as  nodal 
variables.  This  immediately  introduces  the  primary  computational 
disadvantage  of  the  mixed  model  compared  to  the  displacement  model  in 
that  there  are  five  degrees  of  freedom  per  node  for  the  mixed  model 
whereas  only  two  degrees  of  freedom  per  node  in  the  displacement 


model.  However,  this  disadvantage  must  be  weighed  against  the  positive 
aspects  of  the  mixed  model.  The  computed  nodal  stresses  are  typically 
more  accurate  for  the  mixed  model  and  also  appear  as  part  of  the 
solution  vector  instead  of  requiring  any  post  computations  using  the 
gradients  of  the  displacements  and  the  stress-strain  relations  of  the 
material.  This  fact  is  of  immediate  value  in  the  analysis  of  contact 
problems  since  the  boundary  stresses  are  required  to  compute  the  state 
of  stick  or  slip  at  the  nodes  and  also  to  compute  tangential  surface 
tractions  due  to  friction  for  a  contacting  body. 

In  the  purely  numerical  analysis  of  contact  problems,  different 
assumptions  are  made  to  approximate  the  behavior  of  the  contacting 
bodies  at  the  region  of  contact.  For  example,  in  most  analyses,  the 
material  of  one  of  the  bodies  is  allowed  to  penetrate  the  domain  of  the 
second  contact  body  and  is  then  pushed  back  out  during  a  subsequent 
portion  of  the  analysis.  In  addition,  the  coefficients  of  friction  used 
in  the  analysis  are  those  determined  from  physical  tests  of  specimens 
made  of  the  materials  in  contact  and  are  nearly  always  specified  to  be  a 
constant  throughout  the  domain  for  the  duration  of  the  analysis.  These 
are  two  very  basic  but  major  assumptions  frequently  employed  in  the 
development  of  numerical  schemes  used  to  analyze  contact  problems. 
Although  these  assumptions  are  necessary  in  order  to  develop  a 
relatively  efficient  algorithm,  they  may  have  drastic  effects  on  the 
results  of  the  analysis. 

Considering  these  observations,  a  new  method  of  analysis  will  be 


presented  as  a  part  of  the  present  study  that  will  attempt  to  eliminate 
the  effects  of  some  of  the  approximations  made  in  the  analysis  of 
contact  problems.  This  technique  involves  the  combination  of  an 


experimental  technique,  which  determines  the  true  state  of  the 
displacements  on  the  contact  boundary,  with  a  finite  element  technique, 
which  uses  the  Information  provided  from  the  experimental  technique  to 
determine  the  states  of  displacement  and  stress  within  the  rest  of  the 
domain  under  analysis. 

The  major  objectives  of  the  present  study  are  three-fold.  The 
first  objective  Involves  the  development  of  a  mixed  updated  Lagranglan 
formulation  and  corresponding  finite  element  model  of  the  plane 
elasticity  equations  for  large  deformation  analysis.  Second,  the 
formulation  will  then  be  modified  to  develop  several  numerical 
algorithms  for  the  analysis  of  two-dimensional  contact  problems. 

Finally,  the  mixed  finite  element  model  will  be  combined  with  the 
experimental  technique  of  moire  Interferometry  to  form  a  hybrid  method 
of  analysis  for  the  contact  problem  of  a  pin-loaded  plate.  The  goal  of 
this  final  task  Is  quite  different  from  that  of  the  second  task  in  that 
It  is  mainly  being  implemented  to  exclude  the  effects  of  the  aforemen¬ 
tioned  computational  assumptions  on  the  stress  distributions.  On  the 
other  hand,  the  numerical  schemes  are  self  contained  in  that  they  do  not 
depend  on  any  experimental  data,  and  will  attempt  to  model  the  actual 
behavior  of  physical  problems. 

The  following  sections  contain  the  formulation  of  and  examples  for 
the  different  schemes  alluded  to  in  this  section.  Section  2  contains 
the  development  of  a  mixed  updated  Lagrangian  formulation  and 
computational  scheme  for  the  analysis  of  two-dimensional  elasticity 
problems.  This  formulation  is  extended  in  Section  3  to  develop  two 
computational  algorithms  to  analyze  the  contact  problem  of  a  pin-loaded 
plate.  These  approaches  are  developed  separately  to  consider  the 


distinct  cases  of  allowing  for  a  rigid  pin  or  an  elastic  pin.  In 
Section  4,  a  hybrid  experimental/numerical  technique  is  described  that 
combines  the  experimental  technique  of  moire  interferometry  and  the 
numerical  finite  element  method,  again  for  the  analysis  of  the  pin- 
loaded  plate  problem.  The  results  of  several  numerical  examples  from 
the  algorithms  developed  In  Sections  2-4  are  presented  in  Section  5. 
Finally,  summary  and  conclusions  of  the  present  study  are  given  in 
Section  6. 


2.  GOVERNING  EQUATIONS 


2.1  Introduction 

In  this  Section  we  begin  with  the  statement  of  virtual  work  for  an 
arbitrary  solid  body  under  external  load  and  derive  the  variational 
statement  that  will  be  most  convenient  in  applying  the  mixed  finite 
element  model  to  the  problem.  The  finite  element  approximations  for  the 
displacements  and  stresses  are  then  Introduced  Into  the  variational 
statement  for  two-dimensional  bodies  resulting  In  the  final  matrix  form 
of  the  finite  element  equations.  The  required  specification  of  the 
boundary  conditions  Is  discussed  using  the  variational  form  of  the 
linear  equilibrium  and  stress-displacement  relations  of  the  problem,  and 
an  alternative  formulation  of  the  linear  finite  element  matrices  is 
given  as  well.  A  brief  discussion  of  the  required  order  of  polynomial 
approximation  for  the  mixed  finite  elements  Is  also  given. 


2.2  Mixed  Virtual  Work  Formulation 

We  begin  our  derivation  by  first  assuming  that  the  current 
configuration  Is  known  at  time  t  as  well  as  all  equilibrium 
configurations  previous  to  this  time.  We  desire  the  solution,  i.e.  the 
displacements,  strains,  and  stresses,  at  the  configuration  C£  at  time 
t  +  At  (see  Figure  2.1).  Using  the  principle  of  virtual  displacements 
(see  [51]),  we  can  write 


where 


I  Sjs<2elj>dV  *  5<Jf>  ‘  0 

v2 

«(2F)  =  f  2f,su.dV  +  J  2t.su.ds 

V,  1  1  S„  1  1 


(2.1a) 

(2.1b) 
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ontinuumin  Cartesian  coordinates 


the  Cartesian  components  of  the  Cauchy  stress  tensor  In  the 


configuration  C2  occupying  the  volume  V2, 
u.j  =  the  Cartesian  components  of  the  displacement  vector  In 
going  from  configuration  C^  to  configuration  C2, 

^e^.  *  the  Cartesian  components  of  the  infinitesimal  strain  tensor 
associated  with  u^  which  Is  defined  as 


2e(j 


(2.2) 


=  the  Cartesian  components  of  a  point  In  configuration  C2, 

2 

fj  =  the  Cartesian  components  of  the  body  force  vector  measured 
in  configuration  C2. 

2 

t^  *  the  Cartesian  components  of  the  surface  stress  vector 
measured  In  configuration  C2. 

Here  6  denotes  the  variational  operator  and  au^  denotes  the  virtual 
displacement  vector  (l.e.  the  variation  of  u^). 

Although  this  statement  is  valid  at  time  t  +  at.  It  is  not 
immediately  useful  since  the  Integrations  are  performed  over  domains 
that  are  not  yet  known.  However,  we  may  transform  this  equation  into  a 
known  configuration  using  appropriate  stress  and  strain  measures  (52], 
To  do  this,  we  first  define  the  second  Piola-Kirchhoff  stress  tensor  as 


2,  °o  3Xi  2 

1  1j  d  3xm  mn  3xn 


(2.3) 


where  are  the  Cartesian  coordinates  of  a  generic  point  In 
configuration  Cp  pQ  Is  the  density  in  C^  and  p  denotes  the  density  In 
C2.  We  have  indicated  In  the  notation  that  the  second  Piola-Kirchhoff 
stress  tensor  Is  measured  in  C2  but  is  referred  to  Cj.  Since  the  Cauchy 
stress  tensor  is  defined  as  force  per  unit  area  of  the  deformed 


configuration.  It  Is  always  measured  In  and  referred  to  the  most  current 
configuration.  Hence,  we  may  write 


2$  _  2 
2*1j 


T1j  5  2T1 j 


(2.4) 


We  next  define  the  components  of  the  Green- Lagrange  strain  tensor 


as 


?Ei  j 
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au_  au. 
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aX^  aXj 


(2.5) 
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Just  as  the  Cauchy  and  second  Plola-KIrchhoff  stress  tensors  are  related 
by  a  kinematic  transformation  (2.3),  the  variation  In  the  Green-Lagrange 
strain  tensor  and  the  variation  in  the  Infinitesimal  strain  tensor 
e^j  are  related  by 


rtr  \  in  II  .  /  *  \ 

^rij'  '  aX1  aXj  M2emn' 


(2.6) 


Clearly,  both  of  these  strain  tensors  use  the  particle  displacements  u^ 
In  going  from  configuration  C^  to  configuration  C£. 

We  next  note  that  the  second  Plola-KIrchhoff  stress  tensor  Is 
energetically  conjugate  to  the  Green-Lagrange  strain  tensor  and  the 
Cauchy  stress  tensor  Is  energetically  conjugate  to  the  Infinitesimal 
strain  tensor.  In  other  words,  if  we  use  the  definitions  in  Eqs.  (2.3) 
and  (2.6)  along  with  the  Identities 


and 


ax.  aX. 

- —  — “■  s  5 

aX*  ax„  im 
j  m 


(2.7) 


pQdV1  =  pdV2  (2.8) 

then  we  may  write  the  expression  for  the  internal  virtual  work  given  in 


(2.9) 


L  lSij  6(lE1j)dV  =  Sj  6(2e  1j)dV 

V1  v2 


Substituting  Eq.  (2.9)  Into  Eq.  (2.1),  we  obtain  our  modified 
statement  of  virtual  work,  which  Is  now  In  terms  of  a  known 
configuration,  and  Is  given  by 

0  a  i*  lS1j  «<iEij>  dV  -  «(lp>  (2.1°) 

V1 

We  have  also  Implied  the  assumption  that  the  applied  loading  Is 
Independent  of  the  deformation  of  the  structure,  and  hence 


«(ZF)  -  «(1F)  (2.11) 

Next,  we  use  the  Incremental  decompositions  of  the  stresses  and  j 

strains  to  write  I 

I 


isij  *  Sj +  isij 
lE1j  =  le1j  +  ln1j 


(2.12) 


where 


jS.jj  *  Incremental  components  of  2nd  Plola-Kirchhoff  stress  tensor 
^e^ j  *  (Incremental)  components  of  the  Infinitesimal  strain  tensor 


ln1j  = 

Substituting 
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2  3X1  3Xj 


Eq.  (2.12)  into  Eq.  (2.10)  gives 


(2.13) 
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o-J.  <Sj  ♦  is«>  “l'u  ♦  <2-14> 


L  istj  6*ieij  +  inij^  dv  +  L  1,ijs(inij>dv 


*  - 1.  \  5<ieu)dV  +  5(1f)  <2-15) 


We  next  linearize  the  equations  by  assuming  that 


lS1j  *  1^1 jrs  ers  *  6lE1j  =  6le1j 


(2.16) 


and  thus  obtain  the  approximate  governing  equation 


L  1C1  jrs  lers  6(le1j)  dV  +  J'y  Sj  6(ln1j)dV 


*  “  /  6(xei j)dV  +  «(lp)  (2.17) 


The  above  linearization  can  be  Interpreted  as  a  representation  of 


the  nonlinear  curve  between  two  consecutive  load  steps  by  a  linear  line 


segment,  and  must  be  solved  Iteratively. 


A  mixed  (or  stationary)  virtual  work  statement  that  treats  the 


displacements  and  stresses  as  independent  variables  can  be  derived  from 


Eq.  (2.17)  as  follows.  First,  we  note  that  Eq.  (2.17)  is  the  first 


variation  of  the  functional 


n  “  J  1  1C1  jrs  lers  le1j  dV  +  L  Sj^ij  +  lnij)dV  '  ^  (2,18) 
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Next,  we  Introduce  the  stresses  as  additional  dependent  variables  by 


treating  the  strain-displacement  relations 


V  i'.  (Vi*,  i'.T,. 


(2.19) 
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as  constraints.  The  Increments  In  the  2nd  Piola-Kirchhoff  stress  tensor 


components  act  as  the  Lagrange  multipliers  (see  [511).  The 
modified  functional  for  the  mixed  formulation  becomes 


,  3U.  3U  . 

1  /  l  . 


nL  =  n  "  lS1jlle1  j  "  2  (3Xj  +  3X*, 


(2.20) 


where  n  is  given  by  Eq.  (2.18). 


We  next  write  the  linearized  expressions  for  the  strain  energy 


density  UQ  and  the  complementary  strain  energy  density  UQ  due  to  the 


Incremental  displacements  as 


Uo  =  2  Cijk*  le1j  leka 


(2.21) 


Jo  *  1  ui ski  ri  j 


(2.22) 


where  are  the  components  of  the  compliance  tensor.  The  strain 


energy  density  and  the  complementary  strain  energy  densities  are  related 


according  to 


*  Uo  =  Uo  "  lS1j  le1 j 


=  2  C1  jka.  lei j  lekt  "  lSij  leij 


(2.23) 


Using  Eqs.  (2.22)  and  (2.23)  in  Eq.  (2.20),  we  obtain 


1  1  i  ★  1 

■L  •  l 1  'u(ieu +  i«ij) ♦ 1  isu(?4  *  -  uo'dv  -  (*•**> 


Imposing  stationarity  on  this  functional  with  respect  to  the 


displacements  and  the  stresses,  we  obtain  the  two  approximate 


equilibrium  equations 


L  \jfi<l’Mj>dV  +  lVUUdV  *  6(  F)  ’  (2.25a) 

V1  V1  V1 

L  u1,j*(lS1j)dV  "  L  °1 jka  lSka6(lS1j)dV  =  0  (2.25b) 

V1  V1 

These  two  equations  are  analogous  to  Eq.  (2.17)  for  the  displacement 
formulation.  Since  these  equations  are  a  linearized  version  of  the  true 
statement  of  equilibrium,  they  must  be  solved  simultaneously  and 
repeatedly  for  a  given  load  until  the  Increments  of  the  displacement  and 
stress  components  are  within  a  preassigned  tolerance.  The  final 
solution  vector,  after  convergence,  will  represent  the  true  state  of 
equilibrium  for  a  given  load,  and  the  Iterations  are  therefore 
frequently  termed  equilibrium  iterations. 


2.3  Finite  Element  Model 

We  next  construct  the  finite  element  model  of  Eqs.  (2.25)  for  the 
two-dimensional  case.  We  begin  by  assuming  Independent  approximations 
of  the  displacements  and  stresses  of  the  form 

n 


ui(x1.x2)  =  ulVxl,x2) 


(2.26) 


isi j(xi*x2)  =  isijMxrV  (2*27* 

i/ 

where  denote  the  value  of  at  the  k-th  node.  Substituting 
these  expressions  into  the  two  equilibrium  equations  in  Eq.  (2.25)  gives 
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Here  n  is  equal  to  the  number  of  nodes  in  the  element  and  represent 
the  Cauchy  stress  components  that  have  been  determined  at  the  last  known 
configuration. 

As  mentioned  In  the  previous  section,  the  equilibrium  equations  are 
only  approximate  since  we  have  linearized  the  true  equations  of 
motion.  For  this  reason,  there  may  be  errors  introduced  Into  the 
computed  solution  at  each  load  level,  particularly  if  the  load  increment 
is  large.  To  correct  for  this,  we  minimize  the  force  imbalance  that 
results  from  the  linearization  process  for  a  given  load  increment.  We 
do  this  by  updating  the  stiffness  matrices  and  force  vectors  to  account 
for  the  change  in  the  nodal  positions  and  the  Cauchy  stresses  during  a 
given  load  step.  The  iterations  are  continued  until  the  force 
imbalance,  represented  by  the  right  hand  side  vector,  is  reduced  to 
below  some  convergence  limit.  For  example,  the  displacement  and  stress 
component  increments  at  the  (1  +  l)st  iteration  for  the  solution  at  time 
t  +  at  are  calculated  using 

<(KL]  +  tK"-|)1{4}1+1  ■  {F}L  -  {F}*|L  (2.30) 

where  the  superscripts  L  and  NL  denote  the  linear  and  nonlinear 
contributions,  respectively,  and  the  stiffness  and  force  terms  have  been 
computed  using  the  displacements  and  stresses  known  from  the  previous 
iteration  i. 


Equation  (2.30)  must  be  solved  repeatedly  until  the  force  imbalance 


is  reduced  to  below  or  within  a  fixed  tolerance.  This  amounts  to 
measuring  the  percentage  of  the  displacement  and  stress  increments 
(measured  using  the  Euclidean  norm  of  the  incremental  solution  vector) 
relative  to  the  Euclidean  norm  of  the  total  solution  vector.  When  the 
Increase  in  the  displacements  and  stresses  has  been  reduced  to  below  a 
very  small  percentage  of  the  total  solution,  the  approximate  state  of 
equilibrium  has  been  obtained  for  the  given  load  step,  and  the  load  may 
then  be  Increased  or  the  analysis  may  be  terminated. 

The  solution  of  Eq.  (2.30)  allows  us  to  compute  the  total 
displacements  according  to  the  equation 

2u.  =  *u^  +  u.  (2.31) 


In  the  mixed  formulation,  there  is  no  need  to  compute  the  Cauchy 
stresses  using  the  Almansi  strains  as  is  required  in  the  displacement 
formulation.  Since  the  increments  of  the  2nd  Piola-Kirchhoff  stress 
tensor  are  computed  as  nodal  variables,  we  simply  use  our  incremental 
decomposition  of  the  stress  given  in  Eq.  (2.12)  as 


1  4.  < 

Tij  r ij 


(2.32) 


To  obtain  the  values  of  the  Cauchy  stresses  within  a  given  element  as 
required  in  the  computations  of  the  nonlinear  stiffness  matrix  and  force 
vector,  we  may  simply  use  nodal  stress  interpolation,  or 

Tij(xrx2)  =  ^  TijVxrx2)  (2.33) 

When  the  increment  in  the  ,S.  .  terms  is  reduced  to  be  within  the 

A  I  J 

required  tolerance,  we  have  by  definition  of  the  second  Piola-Kirchhoff 


stress  tensor 


and  we  have  clearly  obtained  the  desired  configuration  C£. 

2.4  Some  Computational  Aspects 

If  the  element  stiffness  matrices  shown  in  Eq.  (2.28)  are  assembled 
in  the  form  shown,  the  matrix  [K^]  will  be  a  null  matrix  for  the  first 
iteration  of  the  first  load  step  since  the  Cauchy  stresses  are  zero 
everywhere  in  the  domain  for  the  undeformed  configuration  CQ.  The 
global  stiffness  matrix  counterpart  of  the  submatrix  [K^]  will  also  be 
a  null  matrix,  resulting  in  an  indefinite  system  of  equations.  The  use 
of  pivoting  will  eliminate  this  problem  but  will  generally  destroy  the 
bandwidth  of  the  system  of  equations.  Mirza  [53]  has  suggested 
premultiplying  the  left  and  right  hand  sides  of  the  global  finite 
element  equation  by  the  transpose  of  the  global  stiffness  matrix,  as  in 
the  least  squares  technique,  to  yield  a  positive  definite  system.  In 
the  present  study,  the  first  node  of  the  finite  element  mesh  was 
selected  such  that  both  of  its  displacements  were  specified  to  be 
zero.  This  results  in  the  value  1.0  being  placed  in  the  diagonal 
position  for  the  first  two  rows  of  the  global  stiffness  matrix  during 
the  imposition  of  the  boundary  conditions.  As  the  Gauss  elimination  is 
performed  on  each  row  of  this  matrix,  the  zeros  on  the  diagonal 
corresponding  to  the  remaining  displacement  degree's  of  freedom  are 
eliminated.  Hence,  as  long  as  the  nodes  are  numbered  as  described  here, 
a  conventional  banded  solver  may  be  used  to  solve  the  global  system  of 
equations. 


An  important  note  on  the  order  of  polynomial  approximations  used  in 
the  mixed  model  is  in  order.  The  variational  statements  of  the 
governing  differential  equations  for  problems  in  plane  elasticity  given 
ir.  Eqs.  (2.25)  show  that  the  stress  components  appear  undifferentiated 
whereas  the  displacement  components  u  and  v  are  each  differentiated  once 
with  respect  to  the  x  and  y  coordinate  directions.  Hence  to  ensure 
continuity,  the  stress  components  must  be  approximated  by  at  least  a 
constant  within  the  element  and  the  displacement  components  must  be  at 
least  linear  in  both  x  and  y  within  an  element.  The  approximation  order 
should  also  be  such  that  the  mathematical  definitions  of  the 
displacements  and  stresses  are  accounted  for  and  are  consistent  with  one 
another,  i.e.  the  stresses  are  a  function  of  the  gradients  of  the 
displacements. 

The  continuity  requirements  of  the  variational  statement  of  the 
problem  alone  are  not  sufficient  in  guaranteeing  that  the  mixed  finite 
element  matrices  will  be  acceptable  for  the  analysis  of  a  given 
problem.  Unless  the  displacement  and  stress  approximations  are  of  a 
given  order,  the  element  matrices  will  contain  more  than  the  allowable 
three  zero  eigenvalues  corresponding  to  the  three  rigid  body  modes  for 
two-dimensional  bodies.  This  point  was  examined  in  detail  by  Mirza  and 
Olsen  [20]  who  proposed  and  verified  a  completeness  criterion  that 
restricts  the  choice  of  the  order  of  approximation  for  the  displacements 
and  stresses.  The  completeness  criterion  was  given  as: 

The  strains  from  the  stress  approximations 
should  possess  at  least  all  the  strain  modes 
that  are  present  in  the  strains  derived  from  the 
displacement  approximations . 


i 


When  this  criterion  Is  violated,  the  global  stiffness  matrix  In  the 
mixed  model  will  be  singular  even  after  the  imposition  of  the  boundary 
conditions.  Isoparametric  rectangular  elements  are  used  for  all  of  the 
examples  considered  herein  and,  to  meet  the  requirements  of  the 
completeness  criterion,  only  linear-linear  or  quadratic-quadratic 
combinations  are  used  to  approximate  the  distributions  of  displacements 
and  stresses  In  the  mixed  model. 

A  comment  concerning  the  assembly  of  the  element  equations  In  a 
mixed  formulation  deserves  attention.  The  assembly  of  the  element 
stiffness  matrices  over  the  complete  domain  of  the  problem  to  construct 
the  global  stiffness  matrix  requires  that  the  stresses  be  continuous 
across  each  interface  between  all  elements.  Although  this  assumption  is 
valid  for  many  problems,  there  are  cases  where  the  material  properties, 
and  hence  one  or  more  of  the  three  stress  components,  are 
discontinuous.  Such  an  example  can  be  found  in  the  bending  of  a 
composite  beam.  Clearly,  mixed  elements  should  not  be  used  In  the 
analysis  of  such  problems,  since  the  stresses  will  be  erroneous  at  the 
points  of  material  discontiuity. 

One  way  to  circumvent  this  difficulty  is  to  condense  out  the  stress 
degrees  of  freedom  at  the  element  level  so  that  the  continuity  of  the 
stresses  across  the  element  Interface  is  no  longer  enforced.  To  do  this 
we  recall  from  Eq.  (2.28)  that  the  finite  element  equations  can  be 
written  in  partitioned  form  as 


(KU1  [K12]  l{u})  ({F}) 

(K121T  IK22]  ({S})  ( {0}  j 


(2.35) 


The  second  of  these  two  matrix  equations  can  be  written  as 


(2.37) 


Hence  the  nodal  stress  vector  for  a  given  element  becomes 

{S}  •  -  tK22I_1[K12lT{u} 

Substituting  Eq.  (2.37)  Into  Eq.  (2.36),  we  obtain 

[K11] {u}  -  [K12HK2W2r{u}  =  {F}  (2.38) 

or 

tK*]{u}  *  {F}  (2.39) 

where 

IK*]  *  [K11!  -  IK12] IK22]-1(K121T  (2.40) 

Equation  (2.40)  Is  assembled  as  usual,  and  Is  solved  for  the  nodal 
displacements  (after  applying  the  boundary  conditions).  The  nodal 
stresses  that  were  condensed  out  are  then  computed  at  the  element  level 
using  Eq.  (2.37).  Since  the  stresses  are  no  longer  nodal  variables, 
they  will  be  discontinuous  between  elements. 

Two  points  regarding  the  condensation  procedure  using  linear 
elements  merit  some  discussion.  First,  it  can  be  shown  that  the  matrix 
[K*]  in  Eq.  (2.40)  is  precisely  the  element  stiffness  matrix  derived 
from  a  displacement  formulation.  Second,  computing  the  nodal  stresses 
using  Eq.  (2.37)  yields  exactly  the  same  stresses  as  those  computed 
using  the  procedures  typically  followed  in  a  displacement  formulation, 
l.e.,  computing  the  strains  at  the  node  points  and  then  using  the 
constitutive  relations  to  compute  the  stresses.  As  with  the 
condensation  procedure,  the  stresses  computed  in  a  displacement 
formulation  are  not  continuous  between  elements  due  to  the  discontinuity 
of  the  gradients  of  the  displacements. 

Although  the  condensation  procedure  results  in  the  stiffness  matrix 
derived  from  a  displacement  formulation  for  the  case  of  linear  elements, 


this  equivalence  does  not  generally  hold  for  elements  containing  higher- 
order  approximations.  For  example,  a  quadratic  Isoparametric  element 
will  have  the  exact  same  stiffness  matrices  computed  from  Eq.  (2.59)  as 
from  a  displacement  formulation  only  If  the  element  shape  Is 
rectangular.  If  the  element  sides  are  not  parallel  with  one  another, 
the  entries  In  the  two  stiffness  matrices  will  not  be  identically  the 
same,  although  they  should  be  fairly  close  to  one  another.  As  the 
element  shape  differs  from  that  of  a  rectangle,  the  discrepancies 
increase.  Hence  a  quadratic  element  with  curved  sides  will  possess 
larger  differences  between  the  entries  of  the  stiffness  matrices 
computed  from  the  two  different  approaches  than  will  a  quadratic  element 
with  the  shape  of  a  quadrilateral. 


3.  NUMERICAL  ALGORITHMS 

3.1  Introduction 

In  this  section  we  discuss  two  techniques  for  the  analysis  of  two- 
dimensional  elastic  contact  problems.  Contact  problems  have  a  host  of 
computational  difficulties  since  the  regions  of  contact  are  typically 
not  known  as  a  function  of  any  of  the  parameters  of  the  problem  nor  are 
the  regions  of  relative  stick  and  slip  between  the  two  contacting  bodies 
due  to  the  presence  of  friction.  Most  current  numerical  algorithms  that 
solve  contact  problems  are  relatively  complex  and  use  a  number  of 
Iterative  schemes  to  account  for  the  changing  boundary  conditions  and 
regions  of  contact. 

The  nonlinear  mixed  finite  element  model  described  In  Section  2 
forms  the  cornerstone  of  the  methods  described  in  this  section.  The 
displacement  finite  element  model  has  been  used  almost  exclusively  In 
previous  numerical  analyses  of  contact  problems.  Mixed  elements  provide 
the  immediate  advantage  of  computation  of  stresses  as  nodal  variables, 
which  Is  Ideal  for  contact  problems  since  the  stresses  may  be  obtained 
precisely  on  the  contact  boundary. 

Although  a  number  of  contact  problems  may  be  analyzed  by  the  two 
methods  to  be  described  In  this  section,  the  basic  problem  of  interest 
Involves  a  thin,  pin-loaded  plate  under  a  uniform  in-plane  load,  as 
shown  In  Figure  3.1.  The  plate  may  be  orthotropic  or  isotropic,  and  the 
pin  Is  generally  considered  to  be  isotropic.  The  first  algorithm  allows 
for  the  assumption  of  a  rigid  pin,  and  the  second  algorithm  is  much  more 
general  and  allows  for  an  elastic  pin. 
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Figure  3.1 


Elastic  plate  restrained  by  a  pin  and  subjected 
to  a  uniform,  in-plane  load 
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3.2  Rigid  Pin  Contact  Algorithm 

In  general,  contact  problems  Involve  two  or  more  elastic  bodies 
pressing  against  one  another  under  external  or  internal  load.  In  the 
case  of  a  pin-loaded  plate,  the  bodies  of  Interest  are  the  plate  and  the 
pin.  If  the  assumption  of  a  rigid  pin  is  used,  the  analysis  is 
simplified  considerably.  This  assumption  eliminates  the  need  to  analyze 
the  pin,  which  not  only  provides  a  known  point  of  reference  for  the 
ensuing  contact  (i.e.  the  surface  of  the  pin),  but  also  drastically 
reduces  the  resulting  global  finite  element  system  of  equations  since 
there  Is  no  need  to  discretize  the  domain  of  the  pin.  The  assumption  of 
a  rigid  pin  is  reasonable  if  the  modulus  of  elasticity  of  the  pin  is 
much  higher  than  that  of  the  plate.  Analytical  studies  have  also  shown 
that,  in  the  contact  analysis  of  composite  plates,  pin  elasticity  is  not 
an  Important  variable  and  has  a  relatively  small  effect  on  the  resulting 
stress  distributions  [461 . 

One  simple  and  effective  method  for  analyzing  thin,  orthotropic, 
pin-loaded  plates  was  developed  originally  by  Wilkenson,  Rowlands,  and 
Cook  [47]  and  was  later  refined  by  Rahman  et  al.  [481  to  capitalize  on 
the  computational  advantages  that  arise  from  the  rigid  pin  assumption. 
This  method  uses  three  separate  Iteration  steps  to  account  for  the 
incremental  load  level,  the  contact  process,  and  the  effects  of 
friction.  Both  the  original  and  refined  schemes  use  displacement  finite 
elements.  In  the  load  step  iteration,  the  solution  for  a  given  load 
increment  was  treated  as  a  linear  analysis,  i.e.  the  equations  of  linear 
elasticity  were  used.  The  stresses  in  the  plate  were  computed  for  each 
given  load  level  using  the  original  undeformed  configuration  of  the 
plate  along  with  the  final  displacement  vector. 


The  computation  of  contact  stresses  using  displacement  elements  in 
the  analysis  of  contact  problems  may  create  difficulties  since  the 
contact  action  frequently  results  In  very  large  displacement  gradients 
near  the  region  of  contact.  Since  the  required  stresses  in  a 
displacement  model  are  generally  computed  at  the  element  Interiors  and 
are  then  extrapolated  to  the  contact  boundary,  some  type  of  stress 
smoothing  Is  often  necessary  using,  for  example,  a  local  least  squares 
routine  [54]  or  iterative  improvement  on  the  averaged  nodal  stresses 
[55].  Using  mixed  elements,  this  Is  not  necessary  since  the  stresses 
are  computed  as  nodal  variables  and  no  postcomputation  Is  necessary  to 
modify  the  resulting  nodal  stresses.  It  Is  for  this  reason  that  mixed 
elements  would  appear  to  be  advantageous  over  displacement  elements  for 
contact  problems  since  the  stresses  on  the  boundary  are  required  for 
certain  portions  of  the  analysis. 

The  refined  algorithm  developed  by  Rahman  [48]  uses  a  mixed  polar- 
Cartesian  coordinate  system  to  fix  the  proper  displacements  of  the  nodes 
of  the  plate  that  have  come  in  contact  with  the  circular  pin.  An 
iterative  scheme  is  used  to  ensure  that  all  nodes  that  have  come  in 
contact  remain  in  contact  for  a  given  load  step.  In  other  words,  after 
every  iteration  the  positions  of  all  contact  nodes  of  the  plate  that 
have  previously  come  in  contact  with  the  pin  are  corrected  in  the  radial 
direction  so  that  they  remain  on  the  surface  of  the  pin,  which  is 
actually  defined  as  a  set  of  imaginary  points  that  specify  the  region  of 
no  penetration.  If  the  resulting  shearing  stress  for  a  given  contact 
node  is  larger  than  the  induced  radial  stress  multiplied  by  the  nodal 
coefficient  of  friction,  the  node  is  considered  to  be  sliding,  and  it 
may  subsequently  move  in  the  tangential  direction  of  the  pin. 


Otherwise,  the  node  Is  considered  to  be  sticking  to  the  pin  due  to 
friction,  and  It  Is  fixed  to  an  interpolated  position  on  the  pin  for  the 
remainder  of  the  analysis.  This  iterative  procedure  is  repeated  until 
the  sum  of  the  load  steps  has  reached  the  required  load  level.  The 
details  of  this  method  are  more  completely  described  in  reference  [48]. 

3.3  Elastic  Pin  Contact  Algorithm 

The  assumption  of  a  rigid  pin,  which  is  reasonable  for  cases  when 
the  two  bodies  in  contact  have  a  very  large  difference  in  modulus  of 
elasticity,  is  not  usually  valid  for  contact  between  two  generic 
bodies.  The  algorithm  described  in  the  previous  section,  though  useful 
for  certain  problems,  was  mainly  developed  to  demonstrate  the  use  and 
accuracy  of  the  mixed  finite  element  method  for  the  analysis  of  contact 
problems.  For  general  problems  involving  arbitrary  bodies,  it  is 
necessary  to  revise  the  analysis  to  include  the  effects  of  pin 
elasticity,  which  may  be  significant  if  the  two  bodies  are  of  similar 
constitution.  A  second  computational  algorithm  is  described  below  for 
this  more  general  type  of  contact  problem. 


3.3.1  General  Concepts 

To  account  for  the  complications  arising  from  contact  and  the 
presence  of  friction  between  two  elastic  bodies,  we  add  a  Lagrange 
multiplier  contribution  to  our  original  expression  in  Eq.  (2.2 8)  which 
will  represent  the  summation  of  the  total  potential  of  each  of  the 
contact  forces  acting  at  the  nodes  on  the  discretized  contact 
boundary.  In  addition,  the  kineratics  of  the  elements  of  the  two  bodies 
at  the  contact  Interface  must  be  monitored  such  that  the  nodal 
displacements  are  compatible,  i.e.  the  bodies  must  not  overlap.  We 


therefore  will  eventually  Invoke  stationarlty  of  the  modified  functional 

k 

nr  =  n.  -  l  W,  (3.1) 

L  L  1=1  1 

where  k  represents  the  number  of  the  contactor  nodes  on  the  boundary  and 
W  represents  the  total  potential  for  a  given  contact  force  acting  at  a 
given  contact  node.  This  Idea  was  originated  and  developed  using  a 
displacement  formulation  by  Bathe  and  Chaudhary  [36]. 

To  determine  the  total  potential  for  a  contact  force  at  a  given 
contact  node,  we  consider  the  local  geometry  of  a  contactor  node  K  that 
will  penetrate  the  target  domain  whose  boundary  Is  defined  by  the  nodes 
A  and  B  as  shown  In  Figure  3.2.  In  our  discussion,  we  will  assume  that 
the  two  bodies  have  been  discretized  using  linear  elements.  Although 
the  algorithm  could  be  developed  using  higher-order  elements,  this 
somewhat  complicates  the  analysis,  and  linear  elements  were  used  in  this 
study  because  of  their  relative  simplicity  In  many  phases  of  the 
subsequent  analysis.  For  our  problem,  we  define  the  pin  to  be  the 
contactor  body  and  the  plate  to  be  the  target  body  to  remove  the 
ambiguity  of  analyzing  two  elastic  bodies.  Only  the  nodes  on  the  pin 
will  be  required  to  remain  on  or  outside  of  the  domain  of  the  plate 
during  the  loading  of  the  plate,  while  the  nodes  of  the  plate  are 
allowed  to  be  within  the  domain  of  the  pin.  This  Is  a  key  assumption 
and  requires  some  care  when  modeling  the  problem  to  ensure  that  all 
contactor  nodes  (l.e.  contact  nodes  on  the  contactor  body)  are 
originally  outside  the  target  body. 

In  the  formulation  that  follows,  we  assume  that  (i  -  1)  iterations 
have  been  performed  In  the  quest  for  the  equllibirum  configuration  C2  at 
time  t  +  at  for  a  given  loading.  Ouring  the  last  Iteration  (i-1),  the 
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displacements  of  the  nodes  K,  B,  and  A  have  been  such  that  the  contactor 
node  K  has  penetrated  the  domain  of  the  plate  a  distance  CK,  where 


CK  =  | CK |  =  ak  (3.2) 

and 

CK  =  AKxt  +  AKy3  (3.3) 

The  distance  CK  represents  the  minimum  distance  from  the  penetrated  node 
K  to  the  surface  segment  of  the  target  body,  which  is  defined  as  the 
line  segment  between  the  two  nodes  A  and  B  of  the  target  body. 

The  intermediate  configuration  defined  by  the  updated  Cartesian 
coordinates  of  the  nodes  after  iteration  (i  -  1)  Is  clearly  not  in  an 
acceptable  state  of  equilibrium  since  the  contactor  body  has  penetrated 
the  target  body  at  node  K.  This  node  must  eventually  lie  exactly  on  the 
boundary  of  the  plate,  i.e.  on  the  target  segment  defined  by  the  nodes  A 
and  B,  and  its  relative  position  b  between  nodes  A  and  B  must  be  a 
function  of  the  coefficient  of  friction  between  the  two  bodies.  In 
order  for  this  to  occur,  the  displacements  of  the  nodes  A,  B,  and  K  must 
be  adjusted  accordingly  during  the  following  iteration  and  a  contact 
force  must  develop  at  the  contactor  node  K  as  a  result  of  the 
elimination  of  the  overlap  distance.  We  denote  the  latest  estimate  to 
this  contact  force  at  node  K  in  configuration  as 


2*  (1-1)  -  ,0-1)?  +  ,  (1-1)? 


(3.4) 


This  contact  force  is  equal  to  zero  after  iteration  (i  -  1)  and  is 
developed  as  the  overlap  distance  is  eliminated  during  iteration  (i). 

Although  the  contact  force  acts  alone  at  node  K,  it  must  be 
balanced  by  equivalent  nodal  forces  at  nodes  A  and  B  of  the  target 
segment.  Imposing  moment  and  force  equilibrium  on  the  discretized 


target  segment  defined  by  the  contactor  node  K  and  the  target  nodes  A 
and  B  and  solving  these  equations  simultaneously  gives  the  expressions 
for  the  target  segment  nodal  forces  as 
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(3.5) 


As  these  nodal  forces  are  generated,  the  displacement  increments  at  the 
nodes  A,  B  and  K,  which  are  given  by  auj^,  &u^  and  au^1^, 
respectively,  must  occur  such  that  the  overlap  distance  CK  is  eliminated 
during  iteration  (i).  Hence,  the  total  potential  of  the  contact  force 
at  node  K  may  be  written  as 


«k  -  (241>}T(iu'i)i  -  (M,,}M1-1)> 


(3.6) 


where  the  first  term  is  due  to  the  contactor  body  and  the  remaining 
three  terms  are  due  to  the  target  body. 

Since  the  contact  force  is  originally  equal  to  zero  and  is 
developed  during  the  elimination  of  the  overlap  distance,  we  may  write 
the  incremental  decomposition  of  the  contact  force  at  node  K  as 

2t<'>  =  (3.7) 


We  note  that  if  the  two  bodies  are  in  sticking  contact,  the  contact 
force  components  may  be  written  as  in  Eq.  (3.4).  However,  if  the  bodies 
in  contact  are  sliding,  the  contact  force  component  increment  in  the 
direction  tangent  to  the  target  segment  is  equal  to  zero,  or 


where  n$  Is  the  unit  vector  acting  normal  to  the  target  segment  as  shown 
in  Figure  3.3.  This  expression  is  true  because  ax$  acts  as  a  workless 
constraint  force  and  only  the  normal  contact  force  component  may  act  in 
eliminating  the  overlap  distance.  In  addition,  the  value  of  s  will 
change  during  the  iteration  due  to  the  relative  slip  between  the  two 
bodies.  This  value  change  is  assumed  to  be  negligible  for  each 
iteration. 

In  order  to  impose  stationarity  on  the  contact  functional  n  ,  we 
need  the  first  variation  of  Wk.  Using  the  Equations  (3.5)-(3.7),  we  may 
write  the  first  variation  of  the  total  potential  of  the  contactor  node  K 
due  to  sticking  contact  as 
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where  the  vector  notation  has  been  used  to  imply  that  the  x  and  y 
components  of  the  terms  given  within  the  braces  are  represented  and  the 
T  superscript  represents  the  transpose  of  the  vector.  We  note  that  only 
the  displacement  increments  and  the  contact  force  increments  are  allowed 
to  vary,  but  the  contact  force  components  and  the  overlap  distance 
components  are  fixed  scalar  quantities. 
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Figure  3.3  Direction  of  developed  contact  force  after 
elimination  of  overlap 


Similarly,  we  may  write  the  components  of  the  target  segment  unit 


normal  vector  as 


Ss  *  nsxT  +  nsy3  (3-10) 

and  using  Eqs.  (3.5)-(3.6)  and  Eq.  (3.8)  we  may  write  the  expression  for 
the  first  variation  of  the  total  potential  for  the  contact  force  at  node 
K  due  to  sliding  contact  as 
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Here  we  have  written  out  the  components  of  the  terms  involving  the 
contact  force  increments  to  aid  in  showing  the  origin  of  the  terms  in 
the  sliding  contact  matrix  and  force  vector.  These  are  developed  in  the 
following  section. 


3.3.2  Finite  Element  Matrices 

Our  ultimate  objective  in  constructing  the  total  potential  of  each 
of  the  contact  forces  is  to  impose  stationarity  on  the  modified 
functional  In  Eq.  (3.1)  with  respect  to  the  displacement,  stress,  and 
contact  force  component  increments.  Hence,  we  must  separate  the 
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coefficients  of  the  variations  In  the  contact  force  increments  and  the 
displacement  Increments  of  nodes  A,  B,  and  K  for  the  expressions  In  Eq. 
(3.9)  and  Eq.  (3.11).  The  resulting  contact  equations  can  then  be 
written  in  matrix  form  as 
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where  the  entries  of  the  [K£  ]  and  [K£  )  matrices  and  the  contact  force 
vector  {Zr£1-^}  for  the  cases  of  sliding  and  sticking  contact  are  given 
in  the  Appendix.  These  components  are  then  added  to  the  existing  finite 
element  matrices  resulting  from  the  stationarity  of  nR  (see  Section  2). 
Performing  this  step  allows  us  to  write  the  final  finite  element  matrix 
representation  of  the  stationary  constraint  imposed  on  the  modified 
functional  expressed  in  Eq.  (3.1): 
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The  nature  of  the  contact  matrices  will  depend  on  the  state  of 
stick  or  slip  between  each  contactor  node  and  its  corresponding  target 
segment,  and  the  contact  matrix  and  contact  force  vector  must  represent 
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this  current  state.  The  entries  in  the  contact  matrices  and  force 
vectors  are  shown  to  be  added  to  the  standard  mixed  finite  element 
stiffness  matrix  of  a  typical  element  in  Eq.  (3.13)  as  a  matter  of 
convenience.  It  should  be  understood  that  the  constraint  equations  were 
derived  for  a  generic  contact  node  and  its  target  segment,  and  the 
contact  matrix  and  contact  force  vector  entries  must  only  be  added  to 
the  existing  values  corresponding  to  the  proper  global  degree's  of 
freedom  of  the  standard  stiffness  matrix  for  these  contactor  and  target 
segment  nodes. 

Once  the  proper  contact  matrix  and  contact  force  vector 
corresponding  to  stick  or  slip  for  each  of  the  contactor  nodes  have  been 
added  to  the  global  finite  element  matrix  and  total  force  vector,  the 
solution  procedure  is  similar  to  that  of  the  standard  nonlinear 
analysis.  The  global  system  of  equations  Is  solved  repeatedly  for  the 
Increments  in  the  displacements,  stresses,  and  contact  forces  until  the 
Euclidean  norm  of  the  Incremental  solution  vector  and/or  the  total  force 
vector  are  within  a  preassigned  tolerance.  During  the  iterations,  the 
entries  in  the  contact  matrices  and  the  contact  force  vectors  are 
updated  to  reflect  the  most  current  state  of  the  geometry  and  loading. 
Once  the  solution  has  converged,  another  load  increment  may  be  applied 
or  the  analysis  may  be  terminated. 

3.3.4  Determination  of  Stick  and  Slip 

An  important  facet  of  the  analysis  is  the  determination  of  the 
state  of  stick  or  slip  between  each  of  the  contact  nodes  and  their 
corresponding  target  regions  after  each  iteration.  This  step  indicates 
whether  the  sliding  or  sticking  contact  matrices  should  be  imposed  for 
the  next  iteration.  In  Reference  [361,  which  uses  a  displacement 


formulation,  the  total  distributed  tractions  along  the  contactor 
elements  are  computed  after  several  Intermediate  steps  and  are  compared 
to  determine  the  state  of  stick  or  slip  for  a  given  element.  Using  the 
mixed  formulation,  this  part  of  the  analysis  Is  simplified  considerably 
since  the  state  of  stress  Is  known  precisely  at  the  nodes  and  hence  may 
be  computed  for  each  contactor  node  rather  than  for  an  element  side. 
Using  the  stress  transformation  equations  along  with  the  relative  angle 
(with  respect  to  the  fixed  Cartesian  reference  frame)  of  the  target  edge 
of  the  target  segment,  the  normal  stress  component  on  and  the  shearing 
stress  component  tnt  may  be  computed  for  each  of  the  contact  nodes. 
Clearly,  If  the  two  bodies  are  In  contact,  the  normal  stress  component 
of  the  contact  node  acting  on  the  target  body  should  be  compressive.  If 
we  designate  the  static  or  dynamic  coefficient  of  friction  as  y,  we  say 
that  If,  for  a  given  contactor  node, 

l*ntl  >  »\an\  (3-14) 

then  the  node  Is  In  sliding  contact,  and  If 

v|onl  -  I Tnt I  (3.15) 

then  the  node  is  in  sticking  contact.  It  should  be  noted  that  these 
expressions  are  the  computational  equations  used  to  determine  the  state 
of  stick  or  slip  for  a  given  node.  Physically,  however,  Eq.  (3.14) 
would  be  an  equality  since  the  node  begins  to  slide  as  soon  as  the 
tangential  stress  just  exceeds  the  frictional  capacity  of  the  node. 
Equation  (3.15)  would  then  be  changed  to  a  strict  inequality.  This 
physical  situation  is  a  very  minor  factor,  however,  since  this 
bifurcation  point  would  rarely  be  realized  computationally. 

Specifying  y  to  be  the  static  or  dynamic  coefficient  of  friction  in 
these  equations  will  depend  on  if  the  node  was  sliding  or  sticking 


during  the  previous  Iteration.  This  comparison  Is  made  for  all  nodes  In 
contact  after  each  Iteration,  and  the  correct  contact  matrix  Is 
Implemented  for  the  next  Iteration.  Since  the  state  of  stress  Is  not 
known  for  the  first  Iteration  after  the  nodes  have  come  In  contact,  the 
state  of  sticking  contact  was  assumed  so  that  a  non-recoverable 
tangential  displacement  was  avoided. 

3.3.5  The  Computation  of  Contact  Forces 

In  the  case  of  sticking  contact,  the  Increments  In  the  contact 
force  vector  components  are  both  nonzero  since  a  force  may  develop  in 
both  the  normal  and  tangential  directions  of  contact  along  the  target 
segment  due  to  the  two  bodies  sticking  together.  In  the  case  of  sliding 
contact,  the  only  nonzero  Incremental  contact  force  component  Is  that  In 
the  direction  normal  to  the  target  segment,  which  is  automatically 
accounted  for  by  the  contact  matrices  given  in  Eq.  (3.13).  However, 
there  Is  also  a  force  component  that  opposes  the  relative  tangential 
motion  of  the  contactor  node  due  to  the  presence  of  friction  acting 
along  the  target  segment  containing  the  contactor  node.  Although  this 
force  component  does  not  exist  for  frictionless  problems,  meaning  that 
the  sliding  contact  matrices  may  be  Imposed  for  all  iterations  with  no 
corrections,  this  Is  not  true  in  general.  This  section  addresses  the 
Issue  of  the  computation  of  these  tangential  forces  that  must  be  applied 
to  the  contactor  nodes  with  nonzero  coefficients  of  friction  that  are  In 
sliding  contact  along  their  corresponding  target  segments  for  a  given 
Iteration. 

According  to  the  criteria  given  in  Section  3.2.4  to  determine  the 


state  of  stick  or  slip  for  a  given  contactor  node,  we  observe  that  there 
are  four  possible  combinations  of  stick  and  slip  for  the  triplet  of 


nodes  neighboring  and  including  the  contactor  node  K.  These  four  cases 
are  shown  In  Figure  3.4.  The  first  and  fourth  cases  are  relatively 
simple  to  handle  computationally.  In  case  1,  all  three  of  the  nodes  are 
in  sticking  contact,  and  there  are  no  additional  force  contributions  due 
to  friction  resistance  since  the  nodes  do  not  slide.  In  case  4,  where 
all  nodes  are  sliding,  the  only  forces  that  need  be  applied  In  the 
tangential  directions  are  those  due  to  friction.  To  compute  these 
forces,  we  need  to  compute  the  tangential  tractions  acting  on  the 
contactor  segment  that  corresponds  to  the  total  frictional  capacity. 

This  step  must  be  performed  for  each  of  the  two  contactor  segments. 

To  do  this,  we  must  first  compute  the  normal  and  tangential  surface 
tractions  using  the  nodal  stress  values  computed  during  the  previous 
iteration,  along  with  the  equations 


t  =  o  n  +  t  .  n. 
n  nn  n  n t  t 


*t  *  Tntnn  +  °ttnt 


(3.16) 

For  segment  1,  which  contains  two  sliding  nodes,  the  tangential  surface 
traction  is  given  by 


★ 


K 

+ 


(3.17) 


where  the  subscripts  t  and  n  indicate  the  tangential  and  normal 
directions,  respectively,  and  is  the  static  or  dynamic  coefficient  of 
friction  at  node  1,  depending  on  if  the  node  was  in  sticking  or  sliding 
contact  during  the  previous  iteration.  To  compute  the  tangential  force 
corresponding  to  these  surface  tractions,  we  use  the  formula  for  the 
equivalent  nodal  force  given  for  plane  elasticity  problems  as 

F.  *  §  tYds  (3.18) 

t1  S  1  1 


where  S  represents  the  boundary  of  the  given  element.  Since  we  are 
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Figure  3.4  Possible  combinations  of  nodal  stick  and  slip 
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working  with  linear  elements,  the  frictional  force  corresponding  to  the 

* 

modified  tangential  traction  tt,  which  In  this  case  is  constant  along 
the  segment,  will  be  divided  evenly  between  the  two  nodes  of  the 
segment.  This  procedure  is  repeated  for  the  second  contactor  segment, 
which  borders  the  contactor  node  K  on  the  opposite  side,  and  the  total 
summed  force  of  these  two  frictional  force  components  is  applied  to  the 
contactor  node  K  as  an  external  load  for  the  next  iteration. 

The  remaining  possible  cases  involve  a  contactor  segment  with  one 
node  sticking  and  one  node  sliding.  If  the  neighboring  segment  is 
completely  slipping  (case  2)  or  sticking  (case  3),  the  contact  forces 
for  this  segment  only  are  computed  using  the  procedures  described 
previously  for  the  first  and  fourth  cases.  Hence  the  remaining  case 
that  we  need  to  concern  ourselves  with  is  the  situation  shown  in  Figure 
3.5.  Clearly,  the  basic  problem  involves  accounting  for  the  transition 
between  the  zones  of  stick  and  slip,  including  the  computation  of  the 
tangential  contact  force  due  to  the  appropriate  tractions  along  the 
segment  between  the  sticking  node  and  the  sliding  node. 

In  order  to  determine  the  force  due  to  friction  along  a  segment 
acting  on  the  sliding  node,  we  first  need  to  determine  how  much  of  the 
element  segment  is  in  sliding  contact  and  how  much  is  in  sticking 
contact.  For  the  sticking  node,  the  frictional  capacity  of  the  node 
exceeds  the  corresponding  shear  stress,  while  this  inequality  is 
reversed  for  the  sliding  node.  Since  linear  elements  are  being  used,  we 
assume  that  the  frictional  capacity  and  the  shearing  stress  both  vary 
linearly  along  the  contactor  segment.  Hence  there  is  a  point  m  whose 
position  is  defined  by  the  parameter  a,  which  denotes  the  point  of 
transition  between  the  zones  of  sticking  and  sliding  contact.  Using  the 
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c)  Modified  tangential  surface  tractions 


Figure  3.5  Surface  tractions  for  element  side  with 
sections  of  stick  and  slip 


where  m  represents  the  transition  point.  The  values  of  the  tangential 
tractions  in  the  sticking  regions  remain  the  same  as  computed  using  the 
nodal  stress  values.  Using  this  equation  and  the  modified  values  given 
in  Figure  3.5,  we  may  write  the  final  distribution  of  the  tangential 
traction  t  as  a  function  of  the  length  along  the  segment,  denoted  by  x, 
as 


t  = 


tp  -  t, 

|  >x  +  h 

<  ‘3 


if  0  <  x  <  ad 
if  ad  <  x  <  d 


(3.21) 


where  the  values  of  tj,  t2,  and  t3  are  given  in  Figure  3.5. 


To  compute  the  equivalent  nodal  forces  corresponding  to  these 

tractions,  we  use  Eq.  (3.21)  and  substitute  in  the  value  of  t  given  by 
★ 

Eq.  (3.18)  for  t^.  The  tangential  traction  is  no  longer  a  constant 
across  the  length  of  the  contactor  segment,  and  the  force  will  no  longer 
be  divided  evenly  among  the  two  nodes  of  the  contactor  segment. 
Performing  the  required  Integrations  yields  the  desired  tangential  force 
components  for  the  contactor  nodes  K  -  1  and  K  as 


FK-1  *  *1°^!  "  f )  +  t2od^  _  f )  +  +  “(§  ~  1)] 


2. 

F  °  d  +  t 

hK  Zl  6  r2 


a^d  .  +  j/1  ad\ 

—  +  t3d(?  -  T) 


(3.22) 


The  component  of  force  FK_^  Is  already  incorporated  into  the  contact 
force  acting  at  node  K  -  1  as  part  of  the  analysis.  The  magnitude  of 
the  force  FK,  however,  must  be  added  to  the  tangential  contact  force  at 
node  K  computed  from  the  tractions  acting  on  segment  2.  This  second 
force  is  equal  to  zero  if  the  node  K  +  1  is  not  in  contact,  and  is 
computed  using  Eq.  (3.17)  if  the  node  K  +  1  is  in  sliding  contact. 

As  In  the  case  of  computing  the  state  of  stick  or  slip  for  a  given 
contact  node,  it  is  evident  that  having  stress  as  a  nodal  variable  is 
quite  advantageous.  The  normal  and  tangential  surface  tractions  may  be 


readily  computed  using  the  nodal  stress  values  instead  of  using  the 
nodal  contact  forces  or  other  techniques  to  compute  these  quantities. 


4.  A  HYBRID  NUMERICAL/EXPERMENTAL  TECHNIQUE 


4.1  Introduction 

The  primary  difficulty  in  the  numerical  analysis  of  contact 
problems  is  the  lack  of  knowledge  regarding  the  region  of  contact  and 
the  state  of  stick  or  slip  between  the  two  bodies.  If  the  regions  of 
contact  and  the  state  of  stick  or  slip  are  known  for  a  given  load  level, 
the  analysis  is  greatly  simplified.  A  hybrid  experimental/numerical 
technique  is  presented  in  this  section  that  combines  the  experimental 
technique  of  moire  interferometry  with  the  numerical  finite  element 
method  to  form  a  hybrid  technique.  The  new  method  uses  the  advantages 
of  each  of  the  two  separate  methods  to  create  an  accurate  and  powerful 
tool  of  analysis. 

In  this  section,  the  basic  concepts  of  moire  interferometry,  the 
experimental  portion  of  the  hybrid  technique,  are  presented.  The 
details  of  the  hybrid  technique  are  discussed,  and  the  geometric  and 
material  properties  of  the  physical  specimens  and  the  finite  element 
models  are  given  in  preparation  of  the  presentation  of  the  numerical 
results  given  in  the  following  section. 

4.2  Moire  Interferometry 

A  search  for  the  origins  of  the  French  word  "moire11  would  lead  to  a 
fabric  known  as  watered  silk,  which  displays  varying  patterns  of  light 
and  dark  bands.  For  this  reason,  the  so-called  moire  effect  occurs 
whenever  two  similar  but  not  identical  arrays  of  lines  or  dots  are 
arranged  such  that  one  array  can  be  viewed  through  the  other  [56]. 

Moire  methods  have  been  used  in  the  field  of  solid  mechanics  for  several 
decades  to  analyze  deformed  bodies.  During  recent  years,  the 
sensitivity  of  these  methods  has  been  improved  dramatically. 


Moire  interferometry  is  an  experimental  technique  which  may  be  used 
in  solid  mechanics  problems  to  measure  in-plane  displacements  of 
deformable  bodies  under  external  action.  The  method  is  unique  in  that 
it  is  not  dependent  upon  the  geometrical  or  material  non linearities  of 
the  specimen  under  analysis.  All  moire  techniques  use  two  bar-and-space 
gratings,  one  as  a  reference  and  one  that  is  attached  to  the  specimen. 

As  the  specimen  is  deformed,  the  specimen  grating  deforms  along  with  it, 
and  a  contour  map  of  moire  fringes  is  formed  due  to  the  deformed 
specimen  grating  contrasting  with  the  undeformed  reference  grating. 

Since  the  displacements  in  a  given  direction  are  directly  proportional 
to  the  fringe  order,  the  corresponding  displacements  may  be  computed  by 
the  analyst  by  counting  the  number  of  fringes  on  the  contour  map. 
Stresses  may  then  be  calculated  using  the  gradients  of  the  displacements 
and  an  assumed  constitutive  model,  though  this  can  become  quite  tedious 
and  is  not  as  accurate  as  the  displacements  themselves. 

Moire  interferometry  provides  the  required  sensitivity  that  is 
needed  to  accurately  measure  the  small  in-plane  displacements  of  the 
two-dimensional  pin-loaded  plate  contact  problem.  In  particular,  this 
method  provides  valuable  information  on  the  state  of  the  displacements 
(not  stresses)  at  the  interface  boundary  between  the  pin  and  the 
plate.  This  eliminates  the  need  for  theoretical  assumption  on  exactly 
what  is  occurring  at  the  contact  boundary  that  is  so  common  in  most 
numerical  simulations  of  this  problem.  As  mentioned  previously,  it  is 
the  complex  state  of  contact,  stick,  and  slip  that  makes  this  problem  so 
difficult  to  model  numerically,  and  the  use  of  an  experimental  technique 
to  determine  the  contact  conditions  greatly  simplifies  the  computational 


effort. 


Only  the  salient  features  of  moire  Interferometry  have  been 
discussed  In  this  section,  and  Interested  readers  are  urged  to  consult 
the  more  detailed  references  by  Post  [57-59]. 

4.3  The  Hybrid  Technique 

The  two  basic  components  of  the  hybrid  technique  are  the  nonlinear 
mixed  finite  element  method  and  moire  interferometry.  By  far,  the  more 
demanding  portion  of  the  hybrid  technique  is  the  experimental  portion  of 
the  analysis.  Although  not  discussed  in  this  study,  the  experimental 
details  of  the  hybrid  technique  will  be  reported  by  Joh  [60].  Despite 
the  effort  required  to  obtain  accurate  experimental  data  and  the 
relatively  large  computational  time  involved  in  the  finite  element 
analysis,  the  steps  involved  in  the  execution  of  the  hybrid  method  are 
quite  simple  in  concept  and  are  outlined  as  follows. 

The  exact  displacements  around  the  hole  boundary  of  a  physical 
specimen  may  be  measured  for  a  sequence  of  increasing  load  steps  using 
moire  interferometry.  These  displacement  increments,  along  with  the 
loads  applied  to  the  plate,  are  input  as  a  sequence  of  nonhomogeneous 
boundary  conditions  for  the  simulated  problem  analyzed  by  the  nonlinear 
mixed  finite  element  model.  These  steps  can  be  represented 
mathematically  by  considering  the  global  finite  element  matrix  of  the 
plate  written  in  partitioned  form  as 

IK11]  [K12]  [K13]  UuM  /  {F1} \ 

[K12]T  [K22]  [K23]  j{u2}(  =  |{F2}(  (4.1) 

IK13)1  [K23]T  [K33]  ({a}  )  ({F3}) 

where  the  vector  {u^}  contains  the  degrees  of  freedom  corresponding  to 
known  zero  displacements  (for  example,  along  lines  of  symmetry  of  the 
plate),  {u2}  contains  the  degrees  of  freedom  corresponding  to  the  nodes 


of  the  plate  that  are  In  contact  with  the  pin,  and  {a}  contains  the 
degrees  of  freedom  corresponding  to  the  remaining  unknown  nodal 
displacements  and  stresses  throughout  the  domain  of  the  plate. 

The  imposition  of  known  boundary  conditions  on  the  global  finite 
element  matrix  is  a  relatively  straightforward  task  and  only  the  major 
points  will  be  described  here.  Further  details  of  this  procedure  may  be 
found  in  any  textbook  on  structural  analysis  or  the  finite  element 
method  (e.g.  see  [11]).  The  specification  of  homogenetus  boundary 
conditions,  or  in  this  case  the  known  zero  displacements,  results  in  the 
modification  of  Eq.  (4.1)  which  may  then  be  written  as 

*11]  [0]  [0]  “I  ({Ul}\  ({0}  \ 

[0]  [K22]  [K23]  |{u2}(  =  |{F2}(  (4.2) 

jo]  [k23]t  [k33]J  ({a})  ({F3}) 

where  [I]  and  [0]  represent  the  identity  matrix  and  the  null  matrix, 
respectively,  and  {0}  is  the  null  vector.  This  step  is  typically 
performed  in  all  analyses  of  the  plate  regardless  of  whether  a  hybrid 
technique  or  a  numerical  technique  is  being  used  in  order  to  remove  the 
rigid  body  motion  of  the  plate  and  to  account  for  the  symmetry  of  the 
structure.  The  primary  difficulty  in  the  execution  of  the  numerical 
contact  algorithms  described  in  section  3  is  the  determination  of  the 
vector  {u2}  for  a  given  load  step  such  that  there  is  no  penetration 
between  the  nodes  of  the  pin  and  the  plate  and  that  the  regions  of  stick 
and  slip  are  accounted  for.  Hence  in  a  strict  computational  contact 
algorithm  the  vector  {u2}  is  solved  for  in  an  iterative  fashion  for  a 
given  load  step.  Each  Iteration  involves  significant  computation  In 
order  to  eventually  obtain  a  solution  to  the  problem. 
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In  the  hybrid  technique,  the  vector  {ug}  Is  given  by  the 
experimental  technique  and  contains  the  true  displacements  of  the  plate 
boundary  as  It  has  come  in  contact  with  the  pin.  Imposing  the  known 
displacement  boundary  conditions  on  the  global  system  of  equations  Is 
slightly  more  complicated  than  the  procedure  described  in  Eq.  (4.2)  for 
the  homogeneous  boundary  conditions  since  the  specified  displacements 
are  no  longer  equal  to  zero  and  hence  have  an  effect  on  the  remaining 
equations  of  the  system.  The  force  vector  must  therefore  be  modified  to 
reflect  this  condition.  If  the  vector  {u^}  represents  the  known 
displacement  of  the  contact  nodes  around  the  boundary  of  the  plate 
determined  by  the  experimental  technique,  then  the  modified  global 
matrix  equation  after  imposition  of  all  boundary  conditions  may  be 
written  as 
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We  note  that  the  remaining  force  vector  has  been  modified  to  account  for 
the  nonzero  specified  displacements.  The  components  of  this  vector  may 
be  written  as 


_  |/33  j 
hi  M  *1j  u2 


(4.4) 


where  the  indlclal  subscripts  and  superscripts  indicate  the  appropriate 
vector  or  matrix  entries.  Solving  the  modified  matrix  equation  in  Eq. 
(4.3)  will  result  in  the  correct  boundary  displacements  as  well  as  the 
corresponding  displacements  and  stresses  throughout  the  rest  of  the 


Since  the  solution  of  the  Incremental  displacement  and  stress 

a 

vector  is  an  Iterative  process,  the  vector  {u2}  is  nonzero  only  for  the 
first  Iteration  of  the  solution  procedure  for  the  given  load  step,  and 
is  specified  to  be  zero  for  the  subsequent  iterations.  Hence  for  a 
given  load  step,  the  prescribed  contact  boundary  displacements  are 
applied  to  the  plate  along  with  the  appropriate  force  vector  resulting 
from  the  uniform,  in-plane  load  that  is  applied  to  the  plate.  The 
nonzero  displacements  throughout  the  rest  of  the  plate  along  with  all 
nodal  stresses  are  then  solved  for  using  Eq.  (2.40)  until  acceptable 
convergence  criteria  have  been  met.  As  in  the  case  of  the  rigid  pin 
algorithm  described  in  Section  3.2,  the  pin  is  completely  eliminated 
from  the  finite  element  analysis  and  therefore  only  the  domain  of  the 
plate  is  discretized  and  analyzed. 

The  only  approximations  involved  in  this  technique  are  those  due  to 
the  physical  and  mathematical  limitations  of  moire  interferometry  and 
the  finite  element  method.  The  numerically  solved  problem  is  simply  a 
special  boundary-value  problem  with  specified  displacements,  and  no 
other  assumptions  are  involved. 

Since,  in  theory,  the  one  half  of  the  plate  acts  as  the  mirror 
image  of  the  other  half  of  the  plate,  only  one  half  of  the  plate  domain 
is  modeled  for  the  finite  element  analysis.  This  will  significantly 
reduce  the  computational  time  involved.  Experimental  data  have  shown 
that  the  displacements  around  the  hole  of  the  plate  are  not  exactly 
symmetric,  even  for  an  isotropic  material.  This  is  most  probably  due  to 
the  limitations  of  creating  a  perfectly  symmetric  specimen  and  applying 
a  perfectly  symmetric  load.  The  displacements  are  therefore  averaged  in 
both  of  the  coordinate  directions  of  the  plate  to  yield  one  pair  of 


displacement  values  for  a  given  point  of  the  contact  boundary  for  the 
modeled  half  of  the  plate. 

4.4  Description  of  the  Plate 

In  this  section  the  physical  dimensions  and  material  properties  are 
given  for  the  plate  and  the  pin  used  in  the  experimental  portion  of  the 
hybrid  technique.  Also  given  are  the  different  finite  element  meshes 
used  in  the  numerical  portion  of  the  analysis  of  the  plate. 

A  diagram  of  the  plate  used  in  the  experiments  performed  to 
determine  the  boundary  displacements  between  the  plate  and  the  pin  Is 
shown  along  with  the  dimensions  used  in  Figure  4.1.  The  plate  and  the 
pin  were  both  constructed  of  7075-T6  aluminum  (E  =  10,400  ksl  and  v  = 
0.33).  The  thickness  of  the  plate  was  taken  to  be  0.061  inches.  The 
restraining  pin,  with  a  radius  of  0.3745  inches,  was  fixed  to  a 
structure  exterior  to  the  plate  and  hence  Its  center  did  not  deform 
except  for  a  very  small  displacement  due  to  the  bending  of  the  pin 
fixture.  The  plate  was  loaded  by  means  of  pins  passing  through  the  two 
0.5  inch  diameter  holes. 

Several  assumptions  were  made  in  the  finite  element  modeling  of  the 

plate  domain  for  use  in  the  hybrid  technique.  The  assumption  of 

symmetry  was  used  along  the  length  of  the  plate  in  order  to  reduce  the 

total  number  of  degrees  of  freedom  of  the  problem,  thereby  reducing  the 

cost  of  the  finite  element  analysis.  In  addition,  St.  Venants  principle 
was  invoked  near  the  loaded  end  of  the  plate  to  eliminate  the  localized 
effects  of  the  pins  loading  the  plate.  Hence  only  the  first  8.5  inches 
of  the  plate  were  discretized,  and  the  applied  load  was  assumed  to  act 
as  a  uniform  in-plane  load  to  the  shortened  end  of  the  plate  (i.e.  the 
left  end  of  the  plate  in  Figure  4.1). 


Several  different  mesh  configurations  were  used  In  the  finite 
element  analysis  In  order  to  determine  the  effect  of  element  size  and 
approximation  order  and  also  to  compare  results  for  a  variety  of  domain 
and  variable  approximations.  The  first  mesh  used  Is  shown  In  Figure  4.2 
and  contains  191  linear  Isoparametric  elements  and  228  nodes, 
representing  1140  total  degrees  of  freedom.  As  indicated  by  the  type  of 
element  used  to  approximate  the  geometry,  both  stresses  and 
displacements  were  assumed  to  vary  linearly  in  the  finite  element 
approximations.  Since  the  Inner  boundary  of  the  plate  is  circular,  the 
use  of  linear  elements  introduces  some  domain  approximation  error  into 
the  solution  for  this  mesh  and  any  other  mesh  constructed  of  linear 
elements. 

The  second  finite  element  mesh  used  to  model  the  plate  was 
constructed  of  quadratic  elements  and  Is  shown  in  Figure  4.3.  The  use 
of  quadratic  elements  minimizes  the  domain  approximation  error  in  the 
analysis  of  this  problem.  This  second  mesh  contains  104  elements  and 
367  nodes,  which  corresponds  to  1835  total  degrees  of  freedom.  Although 
the  quadratic  elements  do  represent  an  improvement  over  linear  elements 
in  terms  of  modeling  the  circular  boundary  of  the  plate,  their  use  does 
tend  to  increase  the  bandwidth  of  the  global  finite  element  matrix,  and 
hence  the  computational  effort. 

The  nodes  of  both  finite  element  meshes  used  in  this  example  were 
renumbered  using  the  Cuthi 11 -McKee  ordering  strategy  [651  to  reduce  the 
size  of  the  bandwidth  and  thereby  decrease  the  cost  of  the  analysis.  As 
mentioned  in  Section  2,  the  nodes  must  still  be  numbered  such  that  the 
first  node  has  specified  (zero  or  nonzero)  displacements.  Several 
different  contact  nodes  were  selected  as  the  initial  node  in  the 
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renumbering  scheme  and  it  was  found  that  selecting  the  node  of  initial 
contact  (i.e.  the  node  at  x  =  7.375  and  y  =  0.0)  gave  the  node  numbering 
with  the  smallest  resulting  corresponding  bandwidth. 

The  displacements  around  the  boundary  of  the  hole  obtained  by  moire 
interferometry  are  eventually  expressed  in  Cartesian  components  as  a 
function  of  the  angular  position  around  the  inside  of  the  hole.  The 
displacements  were  given  every  0.5  degrees  for  each  load  step.  In  order 
to  retain  as  much  of  the  accuracy  of  these  displacements  as  was 
possible,  the  contact  nodes  of  the  plate  were  located  exactly  at  a  point 
where  the  experimental  displacements  were  computed,  i.e.  on  the  degree 
or  half-degree.  Oue  to  the  storage  restrictions  of  the  finite  element 
model,  only  a  finite  number  of  nodes  could  be  specified  on  the  plate 
boundary,  which  necessitated  ignoring  many  of  the  data  points  from  the 
experimental  analysis.  Hence  for  the  half-plate  model,  the  contact 
nodes  of  the  linear  element  mesh  were  located  at  the  following  (degree) 
locations:  0.0,  1.0,  2.0,  4.0,  6.0,  8.0,  10.0,  12.5,  15.0,  20.0,  25.0, 
30.0,  35.0,  40.0,  45.0,  54.0,  63.0,  72.0,  90.0,  99.0,  108.0,  117.0, 
126.0,  135.0,  144.0,  153.0,  162.0,  171.0,  and  180.0.  In  an  attempt  to 
incorporate  some  of  the  different  data  point  displacements  and  to 
determine  if  this  range  affected  the  resulting  stress  distributions, 
the  nodes  for  the  quadratic  element  mesh  were  placed  at  slightly 
different  locations:  0.0,  1.5,  3.0,  4.5,  6.0,  7.5,  9.0,  10.5,  12.0, 
14.0,  16.0,  19.0,  22.0,  27.0,  32.0,  38.5,  45.0,  50.0,  55.0,  62.5,  70.0, 
80.0,  90.0,  101.0,  112.0,  123.5,  135.0,  246.5,  158.0,  169.0,  and  180.0. 
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5.  NUMERICAL  EXAMPLES  AND  RESULTS 

5.1  Introduction 

In  this  chapter,  the  results  are  presented  for  a  number  of 
numerical  examples  to  demonstrate  the  accuracy  and  efficiency  of  the 
methods  described  in  the  previous  sections.  The  next  three  sections 
contain  results  for  several  elastic  contact  problems,  and  in  particular 
the  pin-loaded  plate  problem,  using  the  rigid  pin  algorithm,  the  elastic 
pin  algorithm,  and  the  hybrid  experimental /numerical  technique. 

It  is  again  emphasized  that  all  of  the  contact  algorithms  contain 
the  nonlinear  mixed  formulation  described  in  Section  2  as  a  foundation, 
and  the  only  procedures  that  vary  in  these  different  computational 
schemes  are  the  assumptions  and  computations  that  account  for  the 
regions  of  contact  and  the  zones  of  stick  and  slip.  Throughout  the 
presentation  of  the  results,  the  effectiveness  of  having  stress  as  a 
nodal  variable  is  demonstrated  and  highlighted. 

5.2  Rigid  Pin  Algorithm  Examples 

The  contact  algorithm  proposed  by  Rahman  in  [48]  was  implemented 
using  the  three  basic  iterations  of  load,  contact,  and  friction  using  a 
geometrically  nonlinear  formulation  along  with  mixed  finite  elements. 
Here  the  results  of  several  example  problems  involving  contact  between 
an  elastic  body  and  a  rigid  pin  are  presented  not  only  to  demonstrate 
the  accuracy  of  the  algorithm  but  also  to  highlight  the  effectiveness  of 
having  stress  as  a  nodal  variable.  The  results  of  the  example  problems 
are  compared  with  available  analytical  and  numerical  solutions. 


P  E  =  21000  psi 


a)  Geometry  and  material  properties  of  cylinder 


b)  Finite  element  mesh  used  for  quarter  cylinder 


Figure  5.  1  Modeling  used  for  infinite  cylinder  resting  on  rigid  plane 
under  uniform  line  load 
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5.2.1  Infinitely  Long  Cylinder  Under  Uniform  Load 

As  a  first  example  we  consider  an  infinitely  long  cylinder  of 
radius  r  =  1  inch  resting  on  a  rigid  plane  and  under  a  uniform  line 
load.  This  problem  was  modeled  using  the  assumptions  of  plane  strain 
and  a  thickness  of  1  inch.  One  quarter  of  the  circular  domain  was  used 
to  model  the  problem  and  was  approximated  by  84  linear  elements.  The 
mesh  is  shown  in  Figure  5.1.  The  load  was  assumed  to  act  at  the  center 
of  the  cylinder  and  was  applied  in  12  increments  with  the  initial 
increments  smaller  than  the  later  increments.  The  problem  was  modeled 
by  assuming  that  the  cylinder  was  in  contact  with  a  rigid  pin  of  very 
large  radius  (R  =  1000  in.)  to  model  the  rigid  plane.  A  tolerance  of 
0.001  (i.e.  0.1  percent)  was  used  for  the  equilibrium  iterations.  The 
modulus  of  elasticity  used  was  21,000  psi  and  the  Poisson's  ratio  used 
was  0.3. 

The  results  of  the  analysis  are  shown  if  Figures  5.2  and  5.3.  The 
numerical  results  are  compared  with  the  Hertz  analytical  solution 
[631.  Figure  5.2  shows  the  contact  pressure  distribution  plotted 
against  the  distance  from  the  original  point  of  contact  and  corresponds 
to  a  total  applied  load  of  56  lbs.  Figure  5.3  shows  the  load  plotted 
against  the  total  contact  area,  the  data  points  of  which  can  only  be 
determined  when  each  successive  node  comes  in  contact  with  the  pin.  The 
results  from  the  contact  algorithm  appear  to  be  quite  good. 

5.2.2  Orthotropic,  Pin-Loaded  Plate 

The  second  example  considers  a  thin,  orthotropic,  pin-loaded  plate 


with  a  hole  of  radius  r  under  a  uniform  in-plane  load,  similar  to  the 
situation  shown  in  Figure  3.1.  The  modeled  plate  is  shown  in  Figure  5.4 
along  with  the  geometrical  and  material  properties  of  the  plate.  Due  to 


symmetry,  one  half  of  the  plate  was  modeled  using  124  linear  elements 
with  156  nodes  for  a  total  of  780  degrees  of  freedom.  The  material 
properties  given  in  the  figure  are  averaged  properties  from  a  number  of 
species  of  wood.  The  pin  was  assumed  to  be  rigid  and  of  radius  R.  The 
plate  was  loaded  to  a  final  load  of  400  pounds  per  inch  of  plate 
thickness,  and  was  applied  In  18  unequal  Increments.  A  constant 
coefficient  of  friction  of  0.7  was  assumed  at  all  points  of  contact 
between  the  plate  and  the  pin  for  all  load  levels  and  is  a  typical  value 
for  wood  on  steel.  No  equilibrium  iterations  were  performed  for  this 
problem. 

Figure  5.5  shows  the  radial  stress  distribution  as  a  function  of 
the  angular  position  around  the  pin  for  the  nodes  that  have  come  in 
contact  at  the  final  load  step.  These  results  are  compared  with  the 
results  obtained  by  Wilkenson  (64]  using  a  finer  mesh  (385  nodes)  and 
quadratic  displacement  elements.  The  comparison  is  very  good. 

5.3  Elastic  Pin  Algorithm  Examples 

The  method  of  analyzing  contact  problems  developed  in  Section  3.3 
is  much  more  general  than  the  algorithm  described  in  Section  3.2  which 
uses  the  assumption  of  a  rigid  pin.  Although  this  assumption  greatly 
simplifies  the  analysis,  it  also  restricts  the  number  of  problems  to 
which  this  technique  may  be  applied,  since  one  of  the  bodies  must  be 
rigid  and  have  a  circular  shape.  This  second  requirement  is  certainly 
valid  for  the  present  study,  as  the  primary  problem  of  interest  is  that 
of  a  pin-loaded  plate,  but  it  should  be  pointed  out  that  the  scheme 
developed  in  Section  3.3  is  much  more  versatile  and  may  be  applied  to  a 
much  larger  variety  of  contact  problems. 


Radial  stress  (psi) 


Angular  distance,  theta  (deg.) 

Figure  5.5  Radial  stress  distribution  around  hole  of  orthotropic  plate 

using  rigid  pin  algorithm 


5.3.1  Infinitely  Long  Cylinder  Under  Uniform  Load 


As  a  first  example  of  the  so-called  elastic  pin  algorithm,  we 
consider  the  Hertz  contact  problem  analyzed  by  the  rigid  pin  algorithm 
in  the  previous  section.  This  problem  was  modeled  using  the  three 
different  mesh  configurations  shown  In  Figure  5.6.  The  geometry  and 
material  properties  of  the  cylinder  are  identical  to  those  shown  in 
Figure  5.1.  The  rigid  plane  in  this  case  Is  modeled  as  a  square  block 
of  very  high  stiffness  with  its  displacements  specified  to  be  zero  and 
would  be  defined  as  the  target  body  according  to  the  terminology 
introduced  in  Section  3.3. 

Instead  of  applying  the  load  to  the  quarter  cylinder  by  means  of  a 
point  load  along  the  vertical  centerline,  as  done  in  the  previous 
section,  we  instead  require  the  upper  horizontal  mid-plane  of  the 
cylinder  to  deform  a  uniform  amount.  Hence,  we  specify  the  vertical 
displacements  of  each  of  the  nodes  along  this  horizontal  mid-plane  to  be 
a  certain  distance  for  each  successive  load  step.  The  total  load  may 
then  be  computed  by  summing  each  of  the  nodal  forces  corresponding  to 
the  specified  vertical  displacements.  For  this  example,  the  load  was 
applied  In  14  non-uniform  displacement  increments  until  a  final 
displacement  of  0.014  inches  was  reached.  Since  this  problem  involves 
friction less  contact,  sliding  matrices  are  always  Imposed  for  a  node  in 
contact,  even  for  the  first  iteration. 

Figures  5. 7-5.9  show  the  contact  pressure  distributions  computed  by 
the  elastic  pin  algorithm,  shown  by  the  dotted  lines,  between  the 
cylinder  and  the  rigid  plane  for  the  cases  corresponding  to  the  total 
applied  displacements  of  0.006  and  0.014.  We  note  that  the  smaller 
displacement  level  corresponds  to  an  applied  loading  of  approximately  P 


Hertz  solution 


— -» -  Numerical  solution 
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Distance  from  centerline  (in.) 


Figure  5.8  Contact  pressure  distribution  for  infinite  cylinder 
using  elastic  pin  algorithm  and  mesh  2 


>***!Vi 


:  ■  »  M  aw  »«mi  »*  «■«  m  W  H"  SCTC’re’TX'TCTVVl'l  VUUWVVT.V.VVVLVi>^  UV\.>\>li>  ITKVKVHxm 


*  56  lbs.  which  Is  near  the  final  load  level  reached  for  the  same 
problem  analyzed  by  the  rigid  pin  contact  algorithm  and  hence  may  be 
compared  with  the  results  given  In  Figure  5.2.  The  computed  stresses 
are  compared  with  the  Hertz  analytical  solution  which  Is  again 
represented  by  the  solid  line  In  the  figure.  The  agreement  Is  quite 
good  for  both  load  levels  and  Improves  as  the  mesh  is  refined. 

The  stress  distributions  have  been  plotted  on  different  graphs  due 
to  the  change  In  the  Hertz  solution  resulting  from  the  differences  in 
the  total  applied  load.  As  the  mesh  Is  refined,  the  force  required  to 
displace  the  horizontal  mid-plane  of  the  cylinder  decreases  since  the 
cylinder  Is  generally  becoming  more  flexible.  This  change  In  load  is 
relatively  small,  but  makes  a  large  enough  difference  in  the  Hertz 
solutions  to  warrant  separate  figures. 

It  Is  also  of  Interest  to  compare  the  changes  in  the  contact 
pressure  distributions  for  different  frictional  conditions.  The  plot  of 
the  contact  pressure  along  the  rigid  plane  measured  at  the  final  load 
step  for  mesh  3  Is  shown  In  Figure  5.10  as  a  function  of  the  distance 
away  from  the  vertical  centerline  of  the  cylinder  for  various  values  of 
the  coefficient  of  friction,  y.  The  plot  of  the  Hertz  analytical 
solution  is  shown  by  the  solid  line  In  the  figure.  He  note  that  as  the 
coefficient  of  friction  Increases,  the  total  load  required  to  displace 
the  horizontal  centerline  of  the  cylinder  a  uniform  amount  also 
Increases.  Hence  in  Figure  5.10,  even  though  each  of  the  total  applied 
forces  corresponds  to  the  same  uniform  displacement  of  0.014  Inches,  the 
Hertz  solution  is  given  for  the  load  that  is  computed  from  the  finite 
element  solution  for  the  case  of  an  Infinite  coefficient  of  friction. 


Figure  5.10  Contact  pressure  distribution  for  infinite 
cylinder  at  maximum  load  with  different 
frictional  conditions (u  =  coefficient  of  friction) 


The  contact  pressure  distributions  shown  in  Figure  5.10  agree 
fairly  well  with  the  analytical  solution.  As  the  coefficient  of 
friction  is  increased,  the  contact  pressure  also  increases,  especially 
for  the  initial  nodes  of  contact.  The  pressures  computed  from  the 
computational  algorithm  for  all  frictional  conditions  are  generally 
higher  than  the  pressures  computed  from  the  analytical  solution  due  to 
the  stiffening  effect  of  the  cylinder  as  it  deforms.  Although  not  shown 
in  this  figure,  the  contact  area  decreases  as  the  coefficient  of 
friction  increases.  The  contact  area  may  only  be  determined  as  each 
successive  contact  node  on  the  boundary  of  the  cylinder  comes  in  contact 
with  the  rigid  plane.  For  example,  for  the  load  level  corresponding  to 
the  horizontal  centerline  displacement  of  0.011  inches,  the  seventh 
contactor  node  is  in  contact  with  the  rigid  plane  for  the  case  of 
frictionless  contact,  but  this  node  has  not  come  in  contact  for  the 
cases  where  y  =  0.3  and  y  *  ®. 

A  second  means  of  computing  the  total  applied  load  required  to 
displace  the  horizontal  centerline  may  be  Implemented  by  integrating  the 
contact  stresses  along  the  length  of  the  contact  area  of  the  rigid 
block.  This  not  only  provides  a  check  on  the  total  applied  force,  but 
it  also  gives  some  Indication  as  to  the  accuracy  of  the  nodal  stresses 
computed  from  the  mixed  finite  element  model.  The  resulting  forces 
calculated  by  integrating  the  contact  pressures  using  a  consistent 
formulation  for  each  of  the  displacement  increments  are  shown  for  mesh  3 
in  Figure  5.11,  denoted  by  the  dashed  line,  and  are  plotted  against  the 
total  specified  displacement  of  the  cylinder  mid-plane.  These  forces 
are  compared  with  the  forces  computed  from  summing  the  equivalent  nodal 
forces  at  the  points  of  specified  displacement,  shown  by  the  dotted  line 
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Figure  5.11  Total  applied  force  comparison  for  infinite  cylinder 
using  elastic  pin  algorithm 
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In  the  figure.  Also  shown  are  the  forces  computed  from  the  exact 
solution  which  are  denoted  by  the  solid  line.  Clearly,  the  applied 
loads  calculated  from  the  two  methods  within  the  finite  element  solution 
are  in  excellent  agreement  with  the  forces  computed  using  the  pressure 
Integration  being  somewhat  higher  than  those  computed  from  the  nodal 
forces.  Both  of  these  loads  are  larger  than  the  corresponding  exact 
loads  computed  for  a  given  displacement,  and  again  this  is  due  to  the 
stiffening  effect  of  the  cylinder. 


5.3.2  Orthotropic,  Pin-Loaded  Plate 

We  next  consider  the  orthotropic,  pin-loaded  plate  previously 
analyzed  by  the  rigid  pin  algorithm  in  Section  5.2.2.  The  allowance  for 
an  elastic  pin  requires  that  the  domain  of  the  pin  be  analyzed  as  well 
as  the  domain  of  the  plate.  Since  the  elastic  pin  algorithm  was 
developed  using  linear  elements,  this  discretization  introduces  another 
source  of  error  into  the  problem  since  the  approximation  of  the  curved 
boundaries  of  both  the  plate  and  the  pin  must  be  approximated  by  a 
series  of  straight  line  segments.  As  with  any  finite  element 
discretization,  however,  this  approximation  error  decreases  as  the  mesh 
is  refined. 

The  location  of  the  boundary  nodes  at  the  pin/plate  interface, 
though  not  difficult  to  model,  does  require  special  attention  to  ensure 
that  the  limitation  of  straight  line  segment  boundaries  does  not  result 
In  the  nodes  of  the  contactor  body  (i.e.,  the  pin)  being  initially 
within  the  domain  of  the  target  body  (i.e.,  the  plate).  The  mesh 
generation  was  performed  such  that  the  nodes  on  the  boundaries  of  both 
the  plate  and  the  pin  were  located  on  the  correct  geometric  locations  on 


the  appropriate  circular  boundaries,  while  all  nodes  of  the  pin  were 
located  outside  the  domain  of  the  plate. 

A  major  disadvantage  of  this  algorithm  in  the  analysis  of  the  pin- 
loaded  plate  problem  is  that  the  pin  must  be  discretized  as  well  as  the 
plate.  While  this  does  not  increase  the  total  number  of  degrees  of 
freedom  to  an  unreasonable  extent,  it  does  hinder  the  effectiveness  of 
the  bandwidth  reduction  scheme  mentioned  in  Section  4.  As  the  nodes  of 
the  pin  come  in  contact  with  the  target  segments  of  the  plate,  the 
global  influence  of  the  degrees  of  freedom  corresponding  to  the  given 
contactor  node  extends  to  the  corresponding  target  node  degrees  of 
freedom  of  the  plate.  Care  must  therefore  be  taken  in  the  numbering  of 
the  nodes  of  the  pin  and  the  plate  to  ensure  that  the  resulting 
bandwidth  does  not  exceed  the  allowable  maximum  as  new  contactor  nodes 
of  the  pin  come  in  contact  with  the  plate. 

The  domain  of  this  problem  was  modeled  using  three  different  mesh 
configurations  where  the  symmetry  of  the  plate  was  used  to  advantage. 

The  initial  coarse  mesh  is  shown  in  Figure  5.12  along  with  the 
appropriate  specified  zero  displacements.  Most  of  the  subsequent  mesh 
refinement  was  restricted  to  the  region  of  probable  contact  to  increase 
the  number  ofcontact  points  between  the  pin  and  the  plate.  This  region 
Is  shown  by  the  darkened  area  in  Figure  5.12.  The  three  element  meshes 
corresponding  to  this  sub-region  are  shown  in  Figure  5.13.  In  general, 
the  region  of  contact  was  estimated  to  be  within  0°  and  45°,  where  0° 
represents  the  initial  point  of  contact  between  the  pin  and  the  plate. 
The  coarse  mesh  (mesh  1  in  Figure  5.13)  consisted  of  99  elements  and  132 
nodes  and  contained  6  specified  contactor  nodes.  Hence  the  total  number 
of  degrees  of  freedom  corresponding  to  his  mesh  is  672  (132  x  5  +  6  x 
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2).  The  refined  domain  approximations  resulted  in  1085  and  1433  total 
degrees  of  freedom  for  mesh  2  and  mesh  3,  respectively. 

The  plate  was  loaded  to  a  final  load  of  400  pounds  per  inch  of 
plate  thickness,  and  was  applied  In  12  unequal  increments  using  an  error 
tolerance  per  load  step  of  0.01.  The  material  properties  used  for  the 


plate  are  identical  to  those  of  the  orthotropic  plate  shown  in  Figure 
5.4,  but  the  material  properties  for  the  pin  were  taken  to  be  those  of 
steel,  and  were  assumed  to  be  given  by  E  =  29,000  ksi  and  v  =  0.3.  The 
physical  dimensions  of  the  plate  and  the  pin  are  identical  to  those 
given  in  Figure  5.4. 

The  radial  stress  distributions  around  the  hole  of  the  plate  are 
shown  as  a  function  of  angular  position  in  Figure  5.14  for  P  =  0.32  Pmax 
and  in  Figure  5.15  for  P  =  Pmax.  In  Figure  5.15  the  stresses  from  each 
of  the  three  different  meshes  are  represented  by  dashed  lines  and  are 
compared  with  the  results  obtained  by  Wilkenson  [64]  which  are  shown  by 
the  solid  line  and  were  obtained  using  the  assumption  of  a  perfectly 
rigid,  circular  pin.  The  results  are  fairly  reasonable  even  for  the 
coarse  mesh,  although  the  stresses  are  significantly  underestimated  near 
the  initial  points  of  contact.  As  the  mesh  is  refined,  however,  both 
the  shape  of  the  stress  distribution  and  the  maximum  stress  values 
approach  the  solution  obtained  by  Wilkenson,  yielding  fairly  good 
agreement  for  the  final  refined  mesh. 

An  important  point  to  note  is  that  the  maximum  stress  given  by  the 
elastic  pin  algorithm  is  not  given  at  the  point  of  initial  contact.  In 
tact,  the  nodal  stresses  do  not  follow  a  regular  pattern  of  decreasing 
as  the  angular  position  away  from  the  initial  contact  point  increases, 
as  reported  by  Wilkenson.  It  is  unknown  if  the  pattern  shown  is  a 
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Figure  5.15  Radial  contact  stress  distribution  using  elastic 
pin  algorithm  with  P=Pmax 


result  of  the  modeling  of  the  problem,  the  contact  algorithm  used,  or  as 
a  consequence  of  using  mixed  elements.  It  should  be  pointed  out, 
however,  that  the  results  given  in  [64]  were  plotted  using  a  smoothing 
routine,  and  In  the  actual  analysis  of  this  problem,  reported  in  [47], 
this  phenomena  of  increasing  stress  with  increasing  angle  was  recorded 
at  several  points  around  the  boundary  of  the  plate.  Furthermore,  in  one 
of  the  few  other  contact  analyses  that  used  mixed  elements,  this  effect 
is  also  shown  rather  dramatically  in  the  analysis  of  the  Hertz  contact 
problem  [45]. 

In  this  example  problem,  the  results  of  the  elastic  pin  algorithm 
are  compared  with  the  results  obtained  by  means  of  a  different  technique 
using  the  assumption  of  a  perfectly  rigid  pin.  To  simulate  this  rigid 
pin,  the  material  properties  of  steel  were  used.  To  investigate  the 
effect  of  higher  pin  stiffness  on  the  resulting  stress  distributions, 
the  analyses  in  this  section  were  repeated  using  a  pin  modulus  of 
elasticity  one  order  of  magnitude  larger  than  that  of  steel.  Poissons 
ratio  was  kept  as  0.3.  For  both  mesh  configurations,  the  displacements 
near  the  contact  region  decreased  approximately  one  order  of  magnitude 
from  the  displacements  using  the  steel  pin.  This  discrepancy  decreased 
rapidly  away  from  the  contact  region,  with  the  displacements  at  the 
loaded  end  of  the  plate  differing  by  only  one  percent.  More 
Importantly,  the  Increase  in  pin  stiffness  had  a  very  small  effect  on 
the  radial  stress  distributions  as  all  computed  stresses  around  the 
plate  were  within  two  percent  of  one  another.  Hence  within  the  context 
of  this  example,  the  rigid  pin  was  adequately  modeled  by  having  a 
modulus  of  elasticity  approximately  one  order  of  magnitude  higher  than 
that  of  the  plate. 


5.3.3  Aluminum,  Pin-Loaded  Plate 


We  next  consider  a  second  example  problem  Involving  a  pin-loaded 
plate.  The  specimen  used  In  this  example  is  the  aluminum  plate  shown  In 
Figure  4.1.  The  plate  Is  restrained  by  an  aluminum  pin  that  is  assumed 
to  be  of  the  same  thickness  (0.061  Inches)  as  the  plate.  The  mesh  used 
to  model  the  plate  and  the  pin  Is  shown  In  Figure  5.16  and  contains  236 
elements,  286  nodes,  and  has  a  total  of  16  possible  contact  nodes,  which 
accounts  for  1462  total  degrees  of  freedom.  The  mesh  used  to  discretize 
the  plate  is  very  similar  to  the  mesh  of  linear  elements  used  In  the 
Implementation  of  the  hybrid  technique  shown  In  Figure  4.2.  Only  the 
region  of  initial  contact  contains  slightly  fewer  elements  and  nodes. 

In  contrast  to  the  previous  example,  the  initial  clearance  between 
the  pin  and  the  plate  is  very  small  for  this  problem.  Hence  even  for  a 
very  small  initial  load  Increment  several  nodes  will  come  in  contact  at 
the  same  time.  The  algorithm  used  to  analyze  this  problem  was  written 
such  that  the  specified  displacement  Increments  corresponding  to  the 
penetration  of  a  contactor  node  are  imposed  one  at  a  time  to  ensure  that 
the  effects  of  each  of  these  displacement  increments  do  not  subsequently 
alter  the  penetration  distance  of  subsequent  contactor  nodes.  The 
loading  for  this  problem  was  applied  in  6  increments  with  initial  small 
load  steps  giving  way  to  much  larger  load  steps  until  a  final  load  of 
1037  pounds  was  applied  to  the  plate.  The  material  properties  used  for 
this  plate  were  the  same  as  those  given  in  Section  4.4.1.  A  coefficient 
of  friction  of  0.15  was  used  in  this  example. 

Figures  5.17  and  5.18  show  the  radial  and  shearing  stress 
distributions  computed  at  the  end  of  the  third  and  sixth  load  steps, 
respectively.  In  each  of  the  figures,  the  results  are  compared  with  the 
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Figure  5.17  Normal  and  shearing  stress  distributions  for 
aluminum  plate  example  at  3rd  load  step 
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Figure  5.18  Normal  and  shearing  stress  distributions  for  aluminum 
plate  example  at  6th  load  step 


corresponding  results  from  the  hybrid  technique  described  in  Section  4 
(see  Section  5.5).  The  radial  stress  distributions  are  in  quite  good 
agreement  with  one  another  where  the  differences  can  probably  be 
attributed  to  the  linear  domain  approximation  of  the  pin,  the 
approximation  of  the  restraining  pin  using  the  assumptions  of  plane 
stress,  and  the  effects  of  the  relatively  large  load  steps  that  force 
several  nodes  to  come  in  contact  for  a  given  increment  due  to  the  very 
small  clearance  for  this  problem.  The  shear  stress  distribution 
computed  from  the  elastic  pin  algorithm  increases  during  the  initial  15 
degrees  of  contact  and  then  abruptly  drops  to  near  zero.  This  peak  in 
fact  Indicates  the  termination  of  the  zone  of  sticking  contact,  and  all 
nodes  in  contact  past  this  point  are  in  sliding  contact.  The  results 
from  the  hybrid  technique  also  show  this  pattern  to  a  slightly  lesser 
extent,  and  the  shearing  stress  does  not  approach  zero  quite  as 
abruptly. 

We  also  note  that  in  general  the  displacements  in  the  direction  of 
the  load  computed  by  the  elastic  pin  algorithm  are  larger  than  those 
given  by  the  experimental  data  from  the  hybrid  technique.  The 
differences  are  largest  near  the  initial  points  of  contact  and  then 
dissipate  away  from  this  region.  This  behavior  is  almost  certainly 
caused  by  modeling  the  pin  as  a  two-dimensional  body  under  plane  stress 
conditions  when  in  reality  this  fixture  is  much  more  three-dimensional 
than  is  the  plate. 

5.4  Results  of  Hybrid  Numerical/Experimental  Technique 

This  section  contains  the  results  of  the  implementation  of  the 
hybrid  technique  described  In  Section  4.  All  analyses  were  performed 
using  either  the  mesh  of  linear  elements  shown  in  Figure  4.2  or  the  mesh 


of  quadratic  elements  shown  In  Figure  4.3.  The  loading  was  typically 
applied  In  9  increments  where  the  values  used  were  given  as  follows  (all 
values  In  pounds):  225,  467,  667,  855,  1037,  1240,  1458,  1670,  and 
1840.  This  load  was  assumed  to  be  uniformly  distributed  along  the 
bottom  of  the  plate.  All  u  displacements  (i.e.  the  displacements  in  the 
direction  perpendicular  to  the  load)  along  the  centerline  of  the  plate 
were  specified  to  be  zero  due  to  the  symmetry  of  the  plate.  All  nodal 
displacement  increments  around  the  hole  interior  were  specified  using 
the  displacement  values  from  the  moire  analysis,  except  for  the  known 
zero  u  displacements  for  the  two  nodes  located  on  the  centerline  of  the 
plate. 

Several  representative  radial  contact  stress  distributions  are 
shown  in  Figure  5.19.  In  all  plots  containing  the  results  from  the 

hybrid  technique,  the  angular  positions  of  the  data  points  correspond  to 

the  diagram  of  the  truncated  plate  shown  in  this  figure.  The  stresses 
shown  here  correspond  to  the  second  load  level  and  the  ninth  or  final 

load  level.  The  stresses  shown  are  taken  from  the  mesh  of  linear 

elements,  represented  by  the  dotted  line,  and  from  the  mesh  of  quadratic 
elements,  represented  by  the  solid  line.  The  compressive  radial 
stresses  are  at  a  maximum  near  the  initial  nodes  of  contact  near  90°  and 
then  decrease  to  near  zero  for  the  lower  quarter  of  the  hole  boundary. 
This  is  as  expected  since  the  points  of  initial  contact  on  the  plate 
have  normals  that  are  closer  to  the  direction  of  applied  load  than  the 
later  contact  nodes,  and  would  therefore  be  expected  to  balance  the 
majority  of  the  load  in  the  radial  direction.  The  stress  free  condition 
Is  well  shown  in  the  region  from  -90°  to  0°  with  the  small  variations 
most  likely  due  to  numerical  and  experimental  scatter.  In  general,  good 
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Figure  5.19  Radial  contact  stress  distributions  around  hole  boundary 
for  linear  and  quadratic  element  meshes  at  2nd  and 
final  load  steps 


agreement  exists  between  the  stresses  computed  using  the  mesh  of  linear 
elements  and  the  mesh  of  quadratic  elements. 

It  should  be  stressed  that  all  of  the  results  presented  from  the 
hybrid  technique  in  this  section  are  based  on  the  experimental  data  from 
a  single  experiment  and  should  be  viewed  with  this  in  mind.  Since  the 
results  are  so  highly  dependent  on  the  experimental  data,  it  is 
necessary  to  avoid  making  sweeping  generalizations  based  on  the  results 
presented  in  this  section.  Nevertheless,  several  preliminary 
observations  can  be  made  from  these  results.  First,  the  region  of 
contact  is  very  large,  even  for  the  lower  load  steps,  but  never  exceeds 
90  degrees.  This  is  most  likely  due  to  the  very  small  initial  clearnace 
between  the  plate  and  the  pin.  Second,  there  are  obvious  fluctuations 
in  the  compressive  stress  in  the  initial  15  degrees  of  contact  for  both 
levels,  and  in  fact  this  pattern  exists  for  each  load  step.  This  is  a 
consequence  of  an  unexplained  sequence  cf  sign  changes  in  the  v- 
displacement  increments  of  the  experimental  data  that  appear  in  this 
region  for  each  load  step.  These  pecularities  may  be  due  to  surface 
flaws,  a  localized  stick/s  lip  effect,  or  some  other  unexplained  physical 
phenomenon.  Figure  5.20  shows  the  radial  stress  distribution  from  0  to 
90  degrees  for  the  final  load  step  using  the  mesh  of  quadratic  elements 
for  a  series  of  increasing  radii  inside  the  plate  domain.  As  the 
distance  from  the  hole  boundary  is  increased,  the  fluctuations  become 
much  smaller  and  are  nearly  eliminated  altogether  for  a  radius  of  0.6 
inches.  Hence  these  oscillations  are  quite  localized  and  are  most 
likely  due  to  a  surface  effect. 

Figure  5.21  shows  the  shearing  stress  distributions  calculated 
using  the  two  different  mesh  configurations  measured  at  the  final  load 
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Figure  5.21  Shear  contact  stress  distribution  around  hole 

boundary  for  linear  and  quadratic  element  meshes 
at  final  load  step 


step.  The  trends  are  generally  the  same  for  both  meshes  except  for 
several  obvious  kinks  that  are  due  to  the  differences  In  the  data  point 
displacements  used  and  the  differences  in  the  two  meshes.  The  values  of 
the  shearing  stresses  are  quite  small  and  are  In  fact  nearly  an  order  of 
magnitude  smaller  than  the  corresponding  radial  stresses  for  this  same 
load  step.  Although  the  stresses  are  of  the  proper  negative  sign,  they 
remain  negative  and  slightly  nonzero  even  In  the  region  of  supposed  non- 
contact.  It  Is  difficult  to  ascertain  If  these  patterns  are  due  to  the 
actual  behavior  of  the  plate  or  if  either  numerical  or  experimental 
difficulties  are  affecting  the  solution. 

The  circumferential  stress  distributions  for  the  final  load  level 
are  shown  in  Figure  5.22  for  both  mesh  configurations.  These  two  curves 
are  In  very  good  agreement  and,  except  for  the  samll  jumps  In  the  region 
of  Initial  contact  which  are  again  probably  due  to  the  discrepancies  in 
the  v-dlsplacement  data,  this  stress  distribution  is  quite  smooth.  This 
stress  component  is  of  little  interest  In  terms  of  the  contact  problem 
and  Is  mainly  shown  here  for  completeness. 

In  the  preceedlng  discussion,  the  region  of  contact  was  always 
assumed  to  be  determined  by  the  region  of  compressive  stress.  Where  the 
radial  stress  went  to  approximately  zero,  the  contact  region  was  assumed 
to  have  ended.  In  terms  of  the  hybrid  technique,  this  Is  the  only  means 
that  Is  available  for  computing  the  angle  or  region  of  contact  between 
the  pin  and  the  plate  for  a  given  load  step. 

Fortunately,  the  experimental  data  provides  a  second  means  of 
computing  the  contact  area.  A  cursory  glance  at  the  u-displacements 
from  the  moire  results  shows  a  positive  displacement  increment  for  all 
of  the  boundary  nodes  from  90°  down  to  some  angle  <t>.  All  data  points 
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Figure  5.22  Circumferential  stress  around  hole  boundary  at 
the  final  load  step 


past  this  angle  have  a  negative  u-dlsplacement  Increment  up  to  the  final 
boundary  node  located  at  180  degrees.  The  angle  $  decreases  as  the 
loading  Is  Increased.  Intuitively,  one  would  expect  that  the  only 
constraint  that  would  cause  a  positive  u-dlsplacement  for  the  nodes  on 
the  hole  boundary  of  the  plate  Is  the  presence  of  the  pin,  and  that  the 
point  of  the  sign  change  would  correspond  to  the  end  of  the  contact 
region.  This  conclusion  Is  reinforced  by  viewing  the  resulting  radial 
stress  distribution  from  the  hybrid  technique  and  noting  that  the  nodes 
bounding  the  change  in  the  sign  of  the  u-dlsplacement  Increments  also 
bound  the  change  In  the  sign  of  the  radial  stress  component  from 
negative  to  positive  (or  near  zero). 

If  the  above  argument  holds,  one  could  recall  the  boundary 
conditions  of  the  mixed  method  given  in  Eqs.  (2.50)  and  (2.51)  and  note 
that  only  the  displacements  on  the  contact  boundary  need  be  specified 
since  the  remaining  portion  of  the  hole  boundary  Is  stress  free.  Using 
the  criterion  given  In  the  preceedlng  paragraph,  only  those 
displacements  corresponding  to  nodes  thought  to  be  on  the  contact 
boundary  were  specified  using  the  experimental  data.  The  resulting 
radial  and  shearing  stress  distributions  are  shown  in  Figures  5.23  and 
5.24,  respectively  using  the  mesh  of  linear  elements.  The  radial  stress 
distributions  are  in  very  good  agreement  with  one  another,  Indicating 
that  the  change  in  sign  of  the  u-dlsplacement  Increments  most  likely 
does  represent  the  end  of  the  contact  region.  There  is  a  small  but 
sharp  jump  in  the  stress  at  the  node  located  at  0°  which  is  most 
probably  caused  by  specifying  the  displacements  for  one  node  and  not 
specifying  the  displacements  for  the  adjacent  node.  The  numerical 
scatter  is  much  smaller  from  -90°  to  0°  and  the  revised  stress  is 
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Figure  5.23  Radial  contact  stress  around  hole  boundary  with 
varying  regions  of  boundary  displacements 
specified  using  linear  elements 
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Figure  5.24  Shear  contact  stress  distribution  around  hole 
boundary  with  varying  regions  of  boundary 
displacements  specified  using  linear  elements 


90.0 


slightly  more  compressive  from  0°  to  90°,  indicating  that  the  specified 
displacements  somewhat  restrain  the  motion  of  the  plate  using  this 
particular  mesh. 

The  shearing  stress  distributions  are  also  in  good  agreement  except 
near  the  end  of  the  contact  region,  where  the  stress  values  resulting 
from  the  specification  of  only  the  contact  boundary  node  displacements 
dips  slightly  more  negative  and  then  goes  to  near  zero  before  following 
the  original  stress  distribution  for  the  remainder  of  the  hole 
boundary.  Again,  the  shearing  stresses  remain  much  smaller  than  the 
corresponding  radial  stress  components. 

A  final  analysis  was  performed  to  determine  the  effects  of 
geometric  nonlinearity  for  this  particular  specimen.  The  final  load 
step  was  applied  along  with  the  corresponding  specified  boundary 
displacements  as  a  single  load  step  In  contrast  to  the  incremental 
analysis  that  had  been  performed  for  the  preceedlng  examples.  The 
resulting  stresses  and  displacements  were  in  excellent  agreement  with 
the  results  from  the  Incremental  solution,  with  all  values  within  1 
percent  of  each  other  and  most  values  much  closer.  The  single  load  step 
analysis  converged  in  three  Iterations  whereas  the  incremental  analysis 
converged  in  two  iterations  for  each  of  the  9  applied  load  increments. 
Hence  It  appears  that  for  this  example,  the  effects  of  geometric 
nonlinearity  are  quite  small.  This  is  no  doubt  due  to  the  relatively 
stiff  material  being  used  and  the  very  small  initial  clearance  between 
the  plate  and  the  pin. 

In  summary,  the  hybrid  technique  provides  a  useful  alternative  to 
strictly  computational  schemes  for  the  analysis  of  the  pin-plate  problem 
In  that  no  approximations  are  involved  regarding  the  location  of  the 
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contact  region  or  the  regions  of  stick  and  slip.  Since  this  information 


is  provided  by  an  actual  physical  test  there  are  several  unfortunate 


drawbacks  to  this  method  that  have  yet  to  be  resolved.  First,  although 


the  region  of  contact  can  be  computed  using  the  results  of  the  hybrid 


technique  with  what  appears  to  be  good  accuracy,  the  regions  of  stick 


and  slip  are  still  unknown  since  the  relative  displacements  between  the 


points  of  the  pin  and  the  plate  in  the  tangential  direction  are  not 


known.  Second,  the  resulting  radial  and  circumferential  stress 


distributions  on  the  hole  boundary  of  the  plate  are  relatively  smooth 


and  appear  to  be  reasonable.  The  shear  stresses,  however,  are  very 


small  compared  to  the  radial  stresses,  have  no  uniform  pattern  around 


the  hole  boundary  and  are  nonzero  in  the  regions  of  no  contact.  Since 


the  material  used  in  the  experiment  was  isotropic,  this  may  mean  that 


the  extensional  strains  and  eaQ  are  measured  with  good  accuracy  but 

II  0  0 


either  there  is  some  problem  in  experimentally  capturing  the  true 


shearing  strains  or  numerically  modeling  the  shearing  stresses  along  the 


contact  boundary  of  the  plate.  Other  conclusions  must  await  the 


computation  of  the  boundary  stresses  using  the  moire  fringe  patterns  and 


the  results  of  additional  experiments,  which  are  to  be  reported  by  Joh 
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6.  SUMMARY  AND  CONCLUSIONS 

6.1  Summary 

A  mixed  variational  statement  and  corresponding  finite  element 
model  were  developed  for  plane  elastic  bodies  undergoing  large 
deformations  using  the  updated  Lagrangian  formulation.  The  mixed  finite 
element  formulation  allows  independent  approximation  of  the 
displacements  and  stresses  and  both  of  these  quantities  appear  as  nodal 
variables.  This  formulation  was  applied  to  several  plane  elasticity 
contact  problems  to  assess  the  efficiency  and  accuracy  of  this  approach. 

Using  the  mixed  formulation,  two  separate  computational  algorithms 
were  developed  for  the  analysis  of  contact  between  two  bodies 
considering  the  effects  of  friction.  The  first  algorithm  assumed  that 
one  of  the  bodies  was  perfectly  rigid  and  circular  in  shape.  The  second 
algorithm  was  much  more  general  and  allowed  for  two  bodies  of  arbitrary 
shape  and  constitution.  Both  techniques  required  the  use  of  several 
iterative  procedures  to  account  for  the  accurate  computation  of  the 
regions  of  contact  and  the  regions  of  stick  and  slip  betwen  the  two 
bodies.  The  mixed  formulation  yields  stresses  on  the  boundaries  of  the 
two  bodies  in  contact  since  these  components  are  nodal  variables.  This 
facet  is  a  key  advantage  of  this  formulation  since  the  contact  stresses 
are  not  only  part  of  the  desired  solution  to  the  problem  but  are  also 
used  in  several  intermediate  steps  of  the  analysis.  Both  of  these 
algorithms  were  applied  to  several  example  contact  problems,  including 
the  well  known  Hertz  problem  and  several  examples  involving  contact 
between  a  plate  under  in-plane  load  containing  a  hole  and  a  circular  pin 
located  within  this  hole. 
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A  hybrid  numerical /experimental  method  of  contact  analysis  was 


developed  using  the  mixed  finite  element  formulation  and  the 
experimental  technique  of  moire  interferometry.  This  approach  is  unique 
in  that  there  are  no  approximations  required  to  determine  the  contact 
regions  or  the  states  of  stick  and  slip  as  required  by  all  strictly 
computational  algorithms  since  the  displacements  on  the  contact 
interface  are  provided  by  the  moire  analysis.  The  only  approximations 
involved  In  this  technique  are  those  due  to  the  physical  and 
mathematical  limitations  of  the  experimental  and  numerical  techniques. 
This  hybrid  approach  was  applied  to  the  problem  of  an  aluminum  plate 
with  a  hole  restrained  by  a  pin  and  subjected  to  an  in-plane  load. 

6.2  Conclusions 

The  displacements  computed  from  the  mixed  formulation  applied  to 
several  linear  and  nonlinear  plane  elasticity  problems  are  in  good 
agreement  with  values  computed  using  the  more  conventional  displacement 
finite  element  model.  The  stresses  computed  using  the  mixed  method, 
however,  are  not  only  slightly  more  accurate  than  those  computed  from  a 
displacement  formulation,  but  they  are  also  computed  precisely  at  the 
nodes.  This  is  a  valuable  characteristic  when  applying  this  formulation 
to  problems  with  stress  concentrations  or  contact  problems  since  the 
displacement  formulation  typically  requires  that  the  stresses  be 
computed  within  the  elements  which  are  then  somehow  extrapolated  to  the 
boundaries. 

Both  contact  algorithms  yield  reasonably  accurate  results  for  the 
contact  problems  analyzed  in  this  study.  The  stress  distributions  and 
area  of  contact  for  the  Hertz  problem  are  in  very  good  agreement  with 
the  analytical  solution.  The  pin-loaded  plate  examples  are  generally 
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much  more  difficult  to  analyze  in  a  computational  sense.  This  is 
particularly  true  for  the  second  algorithm,  which  must  approximate  the 
domains  of  both  the  plate  and  the  pin.  The  results  of  the  two 
algorithms  are  compared  with  results  reported  by  other  investigators 
using  both  numerical  and  experimental  methods.  The  results  are 
generally  in  good  agreement  with  existing  solutions,  but  at  times  the 
points  of  maximum  stress  computed  by  the  second  algorithm  are  not  at  the 
points  of  initial  contact,  which  is  typically  the  case  for  this  type  of 
problem.  This  might  be  due  to  the  fact  the  approximation  of  the  contact 
boundaries  of  two  circular  bodies  using  linear  elements  necessitates 
careful  modeling  of  these  bodies  to  ensure  that  they  do  not  initially 
overlap.  This  phenomenon  does  not  occur  in  the  analysis  of  the  Hertz 
problem,  where  the  target  body  is  modeled  by  a  straight  line,  and  for 
most  contact  problems  analyzed  here  the  results  were  consistent  with 
those  reported  elsewhere.  Care  must  be  taken  when  applying  these 
contact  algorithms  to  problems  with  known  stress  discontinuities  or 
singularities,  since  by  definition  the  stresses  must  be  continuous  at 
the  nodes. 

The  stress  distributions  computed  from  the  hybrid  algorithm  for  the 
analysis  of  an  aluminum  plate  shows  very  large  regions  of  contact  for 
the  example  problem  considered.  This  is  due  to  the  very  small  initial 
clearnace  between  the  pin  and  the  plate.  There  is  a  very  distinct  jump 
in  the  radial  contact  stress  distributions  near  the  initial  region  of 
contact  for  all  load  steps,  which  appears  to  be  due  to  several 
unexplained  patterns  in  the  displacement  data.  The  shearing  stresses 
are  very  small  relative  to  the  normal  stresses,  and  the  distributions 
have  no  uniform  patterns.  The  contact  stresses  were  computed  by  the 


hybrid  technique  using  two  separate  approaches:  1)  specifying  all 
displacement  Increments  around  the  hole  boundary,  and  2)  specifying  only 
the  displacement  Increments  corresponding  to  the  region  of  assumed 
contact.  The  assumed  contact  region  was  determined  using  the 
displacement  data  of  the  plate  hole  boundary  nodes  In  the  direction 
perpendicular  to  that  of  the  applied  load  direction.  The  results  of 
these  two  approaches  are  In  good  agreement  with  one  another. 

For  the  example  problem  considered,  the  contact  stress 
distributions  are  similar  regardless  of  whether  a  mesh  of  linear  or 
quadratic  elements  Is  used.  Indicating  that  the  mesh  refinement  yields 
little  change  in  the  stresses.  Although  the  contact  stresses  computed 
by  the  hybrid  technique  are  in  fairly  good  agreement  with  the  results  of 
the  second  computational  contact  algorithm,  no  other  comparisons  are 
currently  available. 

Although  the  hybrid  technique  provides  a  unique  way  of  obtaining 
the  stress  distributions  within  the  plate  domain  with  minimal 
computational  effort  compared  with  the  numerical  contact  schemes,  the 
technique  unfortunately  does  not  yield  the  regions  of  stick  and  slip 
between  the  plate  and  the  pin.  Such  information  could  provide  valuable 
information  on  the  bounds  of  the  static  and  dynamic  coefficients  of 
friction  for  the  contacting  bodies.  It  should  be  noted,  however,  that 
the  moire  technique  does  not  Involve  measuring  stresses;  the  stresses 
are  computed  from  the  measured  displacements/stralns  using  a  stress- 
strain  law  that  is  assumed  to  be  valid  for  the  situation. 

The  computational  algorithm  for  the  contact  analysis  between  two 
elastic  bodies  of  arbitrary  shape  could  be  extended  to  higher  order 
elements  to  allow  the  accurate  modeling  of  two  bodies  with  curved 


boundaries  In  contact,  such  as  In  the  case  of  the  pin-loaded  plate. 

This  task  Is  more  complicated  than  It  first  appears  since  general 
expressions  must  be  developed  for  the  contact  force  potentials  and  the 
forces  due  to  friction  considering  curved  element  sides  and  normal 
surface  vectors  that  are  constantly  changing.  This  approach  would  be 
relatively  straightforward  to  apply  to  problems  such  as  that  of  Hertz, 
where  the  target  body  can  be  defined  by  a  straight  line.  This  would 
greatly  simplify  the  formulation. 

The  mixed  finite  element  model  could  also  be  modified  to  Include 
the  effects  of  material  nonlinearity.  This  could  be  of  significance  In 
the  analysis  of  contact  problems  since  large  stresses  commonly  occur  at 
or  near  the  region  of  contact,  and  nonlinear  material  response  could 
have  a  dramatic  effect  on  the  contact  stress  distributions. 

The  hybrid  technique  consisting  of  the  finite  element  method  and 
moire  Interferometry  may  be  applied  to  additional  Isotropic  or 
anisotropic  plates.  It  may  be  of  use  to  consider  specimens  with  a 
larger  Initial  clearance  between  the  pin  and  plate  than  was  used  In  this 
study.  This  would  decrease  the  total  contact  region  and  Increase  the 
region  of  sticking  contact  as  well  as  help  to  minimize  the  effects  of 
minor  surface  flaws  near  the  initial  points  of  contact. 
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APPENDIX 

The  finite  element  equations  (2.28)  can  also  be  written  In  more 
explicit  form  as 
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Note:  The  Cauchy  stresses  are  computed  using  nodal  stress 
Interpolation  and  not  using  the  Almansl  strains.  For  the  first 
iteration,  the  Cauchy  stresses  are  zero  and  hence  (K^l  Is  a  null 
matrix.  This  yields  an  Indefinite  system  and  an  equation  solver  should 
be  selected  accordingly  for  this  first  Iteration. 

The  contact  matrices  and  contact  force  vectors  given  In  Eq.  (3.13) 
can  be  expressed  In  explicit  form  as  follows: 
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The  above  two  vectors  are  the  same  In  the  case  of  sliding  contact. 
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